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Abstract 


Direct numerical simulation (DNS) is used to investigate a time-dependent turbulent wake 
evolving in a stably stratified background. A large initial Froude number is chosen to allow 
the wake to become fully turbulent and axisymmetric before stratification affects the spread- 
ing rate of the mean defect. The uncertainty introduced by the finite sample size associated 
with gathering statistics from a simulation of a time-dependent flow is reduced, compared to 
earlier simulations of this flow. The DNS reveals the buoyancy- induced changes to the turbu- 
lence structure, as well as to the mean-defect history and the terms in the mean- momentum 
and turbulence-kinetic-energy budgets, that characterize the various states of this flow - 
namely the three-dimensional (essentially unstratified), non-equilibrium (or ‘wake-collapse’) 
and quasi-t wo- dimensional (or ‘two-component’) regimes observed elsewhere for wakes em- 
bedded in both weakly and strongly stratified backgrounds. The wake-collapse regime is 
not accompanied by transfer (or ‘reconversion’) of the potential energy of the turbulence 
to the kinetic energy of the turbulence, implying that this is not an essential feature of 
stratified- wake dynamics. The dependence upon Reynolds number of the duration of the 
wake-collapse period is demonstrated, and the effect of the details of the initial/near- field 
conditions of the wake on its subsequent development is examined. 


1 Background and objectives 

This study is concerned with the effects of stable buoyancy upon the evolution of a turbulent 
wake. These effects can be parameterized by the Froude number F = U/NL, where U and 
L are relevant velocity and length scales and N = [—g(dp/dz)/p 00 ] 1 / 2 is the Brunt-Vaisala 
frequency (where g is the gravitational acceleration, z is the vertical coordinate, p is the 
mean background density and p^ is a characteristic/reference density). The Froude number 
is an inverse measure, being small when buoyancy is important and large when it is not. 
For the flow over a wake-generating body, an appropriate Froude number can be defined as 
Foq = 2Uoo/ND , where Uqq and D are the speed and effective cross-sectional diameter of the 
wake-generating body. A local, time-dependent, alternative is the deficit Froude number, 
Fd = 2Ud/NS , where Ud and S are the maximum mean wake- velocity deficit and mean wake 
width. Since the velocity deficit decreases with distance from the body while the wake width 
increases, the wake Froude number strictly decreases with distance. This implies all wakes 
in stratified fluids will be affected by buoyancy at some point in their development. 

Here we consider the weakly stratified case with initially high Froude number (Fd 1), 
such that the wake is largely unaffected by buoyancy for some distance downstream of 
the body. Weakly stratified wakes evolve through three distinct regimes, which Spedding 
[1] referred to as the three-dimensional (3D), non-equilibrium and quasi- two-dimensional 
states. A review of the towed- wake experiments from which this picture derives is given 
by Diamessis, Spedding & Domaradzki [2], while a model that recovers the three regimes, 
along with two others, is presented by Meunier, Diamessis &; Spedding [3]. 

During the 3D regime the high- FA, wake is initially unaffected by buoyancy and thus 
produces fully three-dimensional turbulent motions for a period of time, while the wake- 
defect Froude number remains greater than order 1. The mean and turbulence statistics 
are consistent with unstratified flows during this time and exhibit self-similar behavior 
indistinguishable from that observed in unstratified wakes (Tennekes & Lumley [4] , Bevilaqu 
& Lykoudis [5]). This situation is in contrast to flows with initial Fd of order one or smaller, 
which are immediately altered by stratification, so that the canonical self-similar state is 
never realized. In the extreme initial - Fd -A 0 case, buoyancy suppresses vertical motion 
to the extent that three-dimensional turbulence never appears (Chomaz et al. [6], Riley & 
Lelong [7]). For towed-sphere experiments, the critical value below which the 3D regime 
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does not exist is ^4; cf. Dommermuth et al . , [8]. 

When the initially three-dimensional turbulence in the high-Foo wake begins to ‘feel’ the 
effects of buoyancy, the non-equilibrium (NEQ) (also known as ‘vertical collapse’) regime 
begins. (The term ‘non-equilibrium’ can be viewed as a synonym for ‘intermediate’ or 
‘transitional’, and should not be interpreted literally: it does not imply, for example, that 
in the pre-NEQ state the rates of production and dissipation of turbulence kinetic energy 
balance, since this type of equilibrium does not happen even in the purely unstratified self- 
similar benchmark.) When the wake enters the NEQ phase, it ceases to spread vertically; 
in fact, the height has been observed in some cases to slowly decrease with time during 
this period (Brucker & Sarkar [9]; Diamessis et al. [2]; Lin & Pao [10], the latter for the 
self-propelled case) , although this may be at least in part a consequence of how the height 
is defined. As we will see below, a distinction should be made between the heights based on 
the mean defect (typically measured using a Gaussian curve fit) and on the extent of the 
wake turbulence (revealed for example by shadowgraph or dye visualizations). 

As the wake stops spreading vertically, its lateral growth continues at roughly the same 
rate observed during the 3D phase (Diamessis et al. [2]). (We reserve the term ‘lateral’ to 
indicate the direction aligned with the axis orthogonal to the streamwise- vertical plane.) 
Whether the rate of lateral spreading remains exactly or only approximately constant as 
the flow passes from the 3D to NEQ states is an open question. Brucker & Sarkar [9] cite 
different power- law exponents for the mean-defect width during the different regimes, while 
Diamessis et al. [2] find t 1 / 3 behavior for the entire wake history. 

The NEQ period is also characterized by a large reduction in the rate of decay \dUd/dt\ 
of the maximum mean defect velocity Ud , and also the rate of decay of the turbulence 
kinetic energy. Spedding [1] and others (e.g. Dommermuth et al. [8]; Diamessis et al. [2]) 
have described these as ‘re-stratification’ effects, attributed to the transfer of potential 
to kinetic energy, as the wake turbulence collapses, and heavy/cold fluid parcels fall, and 
light/warm ones rise, towards their original positions. Although this explanation is plausible 
and consistent with the behavior of the mean defect, it is not the only possibility, since the 
mean and turbulence kinetic energy can also be enhanced by alterations to the production 
and dissipation processes (cf. Brucker & Sarkar [9]). We shall present evidence for a picture 
of NEQ dynamics that is not based on the net transfer of potential to kinetic energy, which 
for the turbulence is consistently in the opposite direction and for the mean does not act 
directly upon the streamwise component. 

The NEQ phase is observed to begin at roughly two buoyancy periods (Nt « 2) for a 
wide range of F ^ > 4 (Spedding [1]). Its duration, on the other hand, has been predicted by 
Meunier et al. [3], and demonstrated by the DNS of Brucker & Sarkar [9] and the large-eddy 
simulations (LES) of Diamessis et al. [2], to increase with Reynolds number = U^D /v. 
While the evidence for this Reynolds-number dependence is strong, it remains at present 
rather qualitative, since it is influenced by the statistical uncertainty associated with finite- 
domain numerical simulations (see below). 

After the NEQ phase, the wake enters the quasi-two-dimensional (Q2D) state, char- 
acterized by the emergence of large-scale coherent structures. (Godeford h Cambon [11] 
point out that ‘quasi- two-component’ is a more accurate description of flows like this, whose 
vertical motions have been strongly attenuated by stable buoyancy, since significant vertical 
variations of the other two components are not necessarily negligible. Nevertheless, to be 
consistent with previous studies, we shall continue to use ‘Q2D’ to refer to this aspect of 
the wake evolution.) The towed-sphere measurements of Chomaz et al. [6] revealed that 
flattened ‘pancake-like’ eddies eventually appear for both weak and strong stratification (cf. 
Spedding [12]). The maximum mean velocity begins to decrease more quickly than it did 
during the NEQ regime, while the mean defect spreads laterally at a rate comparable to 
the 3D value. The vertical spreading also begins again, growing at a rate consistent with 
viscous diffusion (Meunier et al. [3], Brucker & Sarkar [9], Diamessis et al. [2]). 
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While the effects of buoyancy that lead from the 3D to the NEQ and then the Q2D 
regimes are well-understood on a qualitative level, there is lesser understanding of the de- 
tailed turbulent processes that drive the regime changes. Our aim in the present study is 
to use numerical simulation to quantify the mechanisms involved, and thereby provide an- 
swers to questions such as the one mentioned above regarding the hypothesized transfer of 
potential to kinetic energy during the NEQ regime. Nearly all wake simulations performed 
to date have made use of the temporal idealization where a fixed streamwise segment of the 
wake is evolved in a reference frame attached to the undisturbed fluid. The weak streamwise 
variation of statistics is ignored so that periodic boundary conditions can be used in the 
streamwise direction. This configuration allows for very efficient numerical procedures and 
has proven to yield results that are in reasonable agreement with laboratory measurements 
(cf. Gourlay et al. [13]). Although we employ the same temporal idealization as those before 
us, we have chosen a longer streamwise domain than has been used in prior studies. Since 
the streamwise direction is assumed to be homogeneous, it forms the basis for computing 
the (time-dependent) turbulence statistics. While relatively short domains require fewer 
mesh points and proportionally less computer time, they provide a smaller averaging base 
at each time; a domain of minimal size may contain only one or two of the large-scale turbu- 
lent structures present at that instant in time. While such a simulation progresses without 
difficulty, there can be large uncertainty, which increases in time, in the computed statistics 
due to inadequate sample. Here we use a streamwise domain size that is between 1.2 and 
5.6 times longer than those employed in previous wake simulations (cf. Gourlay et al. [13], 
Brucker & Sarkar [9], and Diamessis et al. [2]). This extended streamwise domain allows 
for a greater population of ‘large-eddy samples’ and provides less uncertain statistics. This 
strategy leads to resolution of earlier discrepancies, such as that in the lateral spreading rate 
during the NEQ regime reported by [9] and [2]. 

As with the discrepancy in the lateral wake spreading rate, there is no consensus on 
the Reynolds-number dependence of the duration of the NEQ regime. This feature is also 
influenced by statistical uncertainty, and prior results reported by Brucker & Sarkar [9] and 
Diamessis et al. [2] have relatively high uncertainty due to limited large-eddy samples. The 
increased domain size and a separate numerical experiment (discussed in detail below) will 
be used to demonstrate the factors affecting the time at which the NEQ state begins to 
transition to the Q2D regime. 

Finally, the initial-condition dependence of the wake development will be examined, in 
light of the non-universality of the self-similar states observed for late axisymmetric wakes 
in unstratified media. The non-universality of the unstratified case is exemplified by mean 
spreading rates and turbulence-to-mean ratios that vary with the shape of the wake-creating 
object (Bevilaqua & Lykoudis [5]) (although they eventually , at distances corresponding to 
1000s of body sizes downstream, relax to a universal value; see Redford et al. [14], referred 
to hereinafter as RCC). This non-universality has been traced to differences in the near- field 
turbulence structures shed from bodies of different shape. For the stratified counterpart, 
the available data suggest that the structures of the early and late wake are unrelated, with 
pancake eddies ultimately resulting from a wide variety of initial states, some of which are 
non-turbulent or even non- wake-like [2,13,15,16]. Furthermore, Spedding [1] showed that for 
sphere wakes, nominal universality is obtained for the mean-defect history Fd(Q/Eoo when 
weighted by F^ 3 . This scaling was generalized to bodies of various shapes by Meunier & 
Spedding [17] (see also [3]) by accounting, via an effective diameter, for the net momentum 
extracted by the body. The founding assumption of their model will be assessed below. 

The issues raised above will be addressed using DNS of a wake at high and relatively 
high Roq. To put this simulation in context, we mention the earlier numerical studies of 
stably stratified wakes generated by const ant- velocity axisymmetric towed bodies performed, 
for high-Foo, by Gourlay et al. [13] (via DNS), Diamessis et al. [2] (LES) and Brucker & 
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Sarkar [9] (DNS), the latter as part of a study of differences between wakes downstream of 
towed- and self-propelled bodies in high- and low-F^ backgrounds (de Stadler &; Sarkar’s 
[18] DNS includes the effect of body acceleration for moderate stratification). Lower 
cases have been simulated by Dommermuth et al. [8], using LES, and (for a range of 
Prandtl numbers) de Stadler et al. [19], using DNS. Provided it can reach relevant Reynolds 
numbers (and produce accurate statistics), a numerical approach has many advantages, 
supplying detail not available from experiment. For example, Gourlay et a/.’s , [13] DNS 
revealed the vortex-ring structure underlying the Q2D pancake eddies, while the LES of 
Diamessis et al. [2] documented the presence of small-scale Kelvin-Helmholtz instabilities 
during the NEQ regime and its Reynolds-number dependence. Dommermuth et a/.’s [8] 
Iow-Fqo LES quantified the relative importance of the terms in the mean- and turbulence- 
kinetic-energy budgets, during the NEQ and Q2D regimes, that ensue when the turbulence is 
immediately affected by buoyancy. Similar exercises were performed by Brucker & Sarkar [9] 
and deStadler &; Sarkar [18] for their towed- and self-propelled-body flows. 

To summarize, we will use DNS to: 

1. Obtain results with improved ‘eddy sample size’ (via increased streamwise domain 
length), to reduce the uncertainty of the statistics, especially the mean wake spreading 
rate. 

2. Assess more completely how the effects of buoyancy alter the wake’s mean momentum 
and turbulence energy budgets. 

3. Determine whether or not there is a net transfer of potential to kinetic energy during 
the wake collapse in the NEQ regime. 

4. Quantify the Reynolds-number dependence of the duration of the NEQ regime. 

5. Demonstrate how the wake development is affected by the details of the initial con- 
dition used in the numerical simulation, comparing with earlier simulation results to 
explore the extent to which mean-defect histories exhibit universality. 

In § 2, the numerical approach, initialization strategy, and DNS parameters are described. 
Section 3 documents the wake development in terms of flow visualisations, mean- flow statis- 
tics, and the budgets of mean momentum, and mean and turbulence kinetic energy. Closing 
comments and conclusions are summarized in § 4. 

2 Problem formulation, numerical strategy and param- 
eters 

Following Gourlay et al. [13] (cf. [2,8,9,18,19]), we (as in RCC) invoke a Galilean transfor- 
mation, and emulate the actual, spatially developing wake with a parallel- flow, temporally 
developing idealization. The resulting streamwise homogeneity allows efficient DNS with a 
standard triply periodic FFT-based pseudo-spectral method [20]. 

For the high-Froude- number flow studied here, vertical density variations are negligi- 
ble compared to the background reference density p <*>, allowing stable stratification to be 
accounted for with the Boussinesq (divergence- free velocity) approximation. Molecular dif- 
fusion of the (temperature-dependent) density field is matched to that for the velocity field 
by assuming the Prandtl number v jot — 1 (cf. [19]), where v is the kinematic viscosity and 
a is the thermal diffusivity, both of which are assumed constant in space and time (and 
therefore independent of temperature). The pressure variable in the incompressible Navier- 
Stokes equations is eliminated using a velocity- vorticity formulation [21]. A low-storage 
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third-order temporal discretization [22] is applied to the nonlinear terms, while the viscous 
terms are treated analytically [23]. 

The coordinate x — x\ denotes the streamwise/longitudinal direction, along the wake 
axis, with u = u\ the corresponding streamwise velocity, defined such that positive values 
point ‘upstream’ and the maximum occurs at the centerline (figure 9). Because of the 
parallel- flow assumption, the uniform background/freestream velocity amounts to an 
arbitrary Galilean shift; we chose Vqo = 0, which implies u 0 far from the wake centerline. 
The lateral (‘spanwise’) and vertical directions are respectively y — x 2 and z — x 3 (positive 
upwards), with velocities v — U 2 and w = U 3 , such that x-y planes are normal to the gravity 
vector Qi = -gS i3 . 

Before discussing the numerical parameters, we describe the three types of averaging used 
in this paper. The Reynolds decomposition is defined with respect to the mean (denoted by 
an overbar) formed by averaging over the homogeneous streamwise direction x\ this mean 
thus varies with lateral and vertical position, and time: e.g. u = u(y,z,t). (Note that, 
unlike for the unstratified wake simulated in RCC, averaging over the azimuthal direction is 
not appropriate.) A prime (') indicates a deviation from this mean, such that the Reynolds 
stresses, density- fluctuation variance and turbulent buoyancy flux are written respectively 
as —UjUj, p' p' and —p'w f . As in Spedding [1], we also examine histories of bulk averages 
of Reynolds- averaged quantities, taken across the inhomogeneous directions (y and z) of 
the mean defect, over the region where u > 0.05 Ud (Spedding’s mean was over horizontal 
and vertical planes where u > 0.2 Ud)- These averages are indicated with angle brackets: 
(k) = (1/2) (uX) = (l/2)f 0 A ‘f 0 A ^u'Xdydz/V, where = ^{v,z)\u/u d>0 . 0 5 is unity if 

u(y,z) > 0.05 Ud and zero otherwise, and T =# J Q Az J Ay i/j dy dz. A third type of statistic 
(indicated by braces { }) is formed from cross- wake integration (not averaging) over the 
entire y-z domain, of Reynolds-averaged quantities; examples include the total/integrated 
turbulence kinetic energy (TKE) {k} = (1/2) f^ z J^ z u' i u' i dydz and its rate of production 
{Vtk e } = ~ ltfo Z (' Ku'jdui/dxj ) dy dz. 

Based on our earlier experiences with DNS of unstratified turbulent wakes (RCC), the 
wake is initialised as a sequence of vortex rings, identical to Case VR of RCC apart from 
the ring spacing (see their figure la), embedded in a stably stratified background with 
imposed constant/uniform mean-density gradient T = ( dp/dz ) OQ < 0. This strategy was 
found to be an efficient way to produce fully developed wake turbulence containing mature 
turbulent structures, which - while not necessarily corresponding to those generated by a 
realistic wake-creating body in a real-world application - allows us to perform a numerical 
experiment that can be used to address the issues raised above (such as the effect of the 
details of the initial conditions on late- wake development). 

In order to simulate the wake using periodic boundary conditions, the transported density 
variable (= ~p(z, t) + p'(x, y, z, t)) is defined as the perturbation from the background density 
field poo + zY (where the origin of the z - axis is, say, in the center of the computational 
domain) . 

As in Shariff [24], each ring in the sequence is represented as a circular loop, of radius ro 
in y-z planes, of Gaussian azimuthal core vorticity u/#, with ujq = (7/71^) exp [— (s/J w ) 2 ], 
where 7 is the ring’s circulation, s 2 = (pc — xq ) 2 + (r — ro) 2 , r 2 = (y — yo ) 2 + (z — z o) 2 and 
( xo,yo,zo ) is the location of the centre of the ring in question. The length-scale 5 U defines 
the thickness of the vortex core (chosen as 6 U = 0.4ro). The rings are spaced uniformly every 
27rro/3 along the wake axis (x,yo 5 ^o) (RCC used 167rro/27 ~ 0.597rro); the rings combine 
to induce an initial mean velocity defect of maximum velocity Udo and width h Q , the latter 
defined as the radial distance from the wake centreline to the location at which the mean 
streamwise velocity u is half of Ud- 

The wake turbulence is the result of a nonlinear /bypass transition, which is triggered 
by a A r = =L0.004ro maximum radial displacement of the ring circumference, in the y-z 
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plane, for each of the azimuthal modes 5 — 32, with random azimuthal phases, using dif- 
ferent random- number seeds for each ring. As the individual vortex rings break down and 
disturb their neighbours, a vericose-type instability ensues (see figure 5a of RCC), leading 
to fully developed turbulence relatively quickly - much sooner than, for example, when 
very-small-amplitude disturbances are added to a mean defect profile, such that transition 
is via growth of the helical linear-instability modes (as was done for Case SD of RCC). The 
approach used by Gourlay et al [13] also avoids the need to simulate the natural linear- 
instability transition process, by prescribing a Gaussian mean velocity and superimposing 
random finite-amplitude perturbations whose magnitudes were set by matching a desired en- 
ergy spectrum, and then clipped/filtered in physical space with an axisymmetric/cylindrical 
window. Similar procedures were followed by Brucker & Sarkar [9] and Diamessis et al [2], 
although the former’s DNS involved an initial adjustment period during which the mean 
flow was held fixed as the fluctuations developed, and the latter’s LES required a precursor 
‘relaxation’ simulation, during which the effective gravity was slowly increased. Despite 
being necessarily somewhat ‘synthetic’ in nature, those initialization schemes all produced 
statistics that very quickly exhibit expected behavior (see figure 21). The primary difference 
between the present and earlier strategies is that here coherent structures are ‘hard wired’ 
into the initial-condition turbulence. 

The initial Froude number specified here is quite large, to ensure the wake turbulence 
is fully developed and physical when it is first affected by the stratification. The initial 
mean-defect Froude number Fdo = Ud 0 /Nh 0 = 23.7; this corresponds to = 2Voo/ND of 
approximately 130, assuming the uniform flow speed Voo ~ 5 U do and virtual-body diameter 
D « 2 h Q (see RCC). The initial mean-defect Reynolds number Rdo = Udoh 0 /v — 2625, 
such that the equivalent Roo — V^D/v ~ 2.6 x 10 4 . This is well above the critical value 
of Roq w 5 x 10 3 , found by Spedding, Browand & Fincham [25] to be required for the 
turbulent wake structures and decay rates to exhibit self similarity - although it is not 
above the threshold (if one exists) beyond which all features (such as the duration of the 
various self-similar regimes) can be expected to be independent of Reynolds number (cf. 
Diamessis et al. [2]). 

A modified version of this case will also be considered, which is equivalent except that 
beginning from Nt = 1.9 the viscosity varies at each time step, such that the mean-defect 
Reynolds number, based on the integral scale r* = ((ln2/7r )(Id/Ud)) 1 ^ 2 , where Id = {u} 
is the volume- flux deficit, remains fixed thereafter at Udr* / v = 1500; see + symbol in 
figure lc, d. (The regridding/domain-expansion strategy described below was also employed, 
although, towards the end of the run, not at identical times.) Results from this ‘artificial’ 
numerical experiment (cf. [26]), which we shall refer to as Case R1500, will be compared 
to those from the primary (constant-viscosity/variable-Renolds-number) DNS, to determine 
the Reynolds-number dependence of the mean-defect histories; unless otherwise stated (see 
figure 7) all the new simulation results shown below will be from the primary DNS. Two 
additional test runs were also performed, to explore numerical- validation issues, as explained 
below. 

For Gourlay et al. [13], Fdo = 34.0 and Rdo = 1472. The more recent simulations tended 
to use lower F and higher R: Brucker & Sarkar’s [9] Case TR50F04 (i.e. F 0 c /2 = Voo /ND = 
4, Roo = 50 000) corresponds to Fdo ~ 0.75 and Rdo ~ 3250 (they assume U do /Voo = 0.11 
and D ~ 1.7 h Q ). They also simulated a lower-R, higher- F flow with which we will compare 
(their ‘TR10F20’), with ( Fd 0 ,Rdo ) ~ (3.75,650). Comparison will also be made with the 
recent F 0 c = 64 (Fdo « 14.5), R 0 0 = 10 5 (Rdo ~ 14 X 10 4 ) LES of Diamessis et al. [2], for 
which the eddy-resolving capability was comparable to the present DNS. In fact, although 
the smallest resolved scales are not those that actually appear at R 0 0 = 10 5 , the LES has 
the advantage of taking full benefit of the available spatial resolution - since the magnitude 
of the subgrid-scale viscosity of the LES can be smaller than the molecular viscosity of the 
DNS - and thus has the capacity to capture scales smaller than those of an over-resolved 
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Figure 1. Mean (a, b) Froude number Fd = Ud/Nh and (c, d) Reynolds number Rd — 

Udh/v histories: , present DNS [{Fdo,Rdo) = (23.7,2223)]; +, time at which the 

variable- viscosity run (Case R1500) first reaches its minimum Reynolds number, U d r* jv — 
1500 (see discussion regarding figure 7); o, Gourlay et al. [13] [( Fd 0 ,Rdo ) ~ (34, 1472)]; ■, 
Case TR50F04 of Brucker & Sarkar [9] [( F do ,Rdo ) ~ (0.75,3250)]; •, Case TR10F20 of 
Brucker & Sarkar [9] [(F^ 0 , Rd 0 ) ~ (3.75, 650)]; o, Case R100F64 LES of Diamessis et al. [2] 
[(Fdo,Rdo)**( 14.5,11000)]. 
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DNS at smaller Reynolds number, such as is used here. 

The variations of the mean-defect-based Froude and Reynolds numbers for the present 
and earlier simulations are shown in figure 1, defined in terms of both the wake’s mean height 
h z and width h y . Evidence of each of the 3D, NEQ and Q2D regimes can be seen. The 
tendency for all the Froude numbers to fall monotonically indicates that all the flows, even 
those with very large F dcn eventually are dominated by buoyancy effects. The anisotropy 
introduced by the stable stratification can be inferred in the differences between the Froude 
and Reynolds numbers based on h z and h y . The wake evolution is such that the power-law 
decay of the height-based Reynolds number Udh z /v varies little over the life of the flow, 
while U d h y /v remains relatively constant, especially during the NEQ phase (cf. Diamessis 
et al. [2]). The time required (t ~ 1/JV) for the vortex-ring initialization (solid line) to 
produce a fully turbulent state is apparent, which can be contrasted with the approaches 
used to define the Nt = 0 state for the other simulations (recall the difficulties involved in 
specifying finite disturbances typical of the fully developed wake downstream of a sphere). 
The question of appropriate virtual origins, needed when comparing the present and earlier 
cases, will be taken up below. 

To prepare for the possibility of the wake generating significant internal/buoyancy waves, 
a ‘fringe’ region is adopted at the lateral and vertical edges of the periodic domain, to prevent 
the energy radiated in the waves from re-entering the domain, and spuriously affecting the 
wake dynamics. A forcing term F (similar to that used by de Stadler et al. [19] and 
de Stadler & Sarkar [18]) is added to the right-hand-side R of the momentum and scalar 
equations 


<9q 

dt 


R + F, 


(i) 


which, in the fringe zone, drives q = (u,v,w, p') T to a mean reference q, defined as its 
average along the edge of the domain. The fringe forcing is written as F = — cr(y, z )( q — q) 
where the fringe function 


a(y : z) = C a [(7 y (y) + a z (z) - a y (y)a z (z)\ , 


(2) 


for 0 < 2 / < and 0 < z < A z , with C a = 0.04 and 


<Ty{y) = 1 - \ ( 

= 1 ~\{ 


tanh 

tanh 


a f 


r - 1 

Ay 

f 2Z 1 

a/ Ivr 


+ tanh 
+ tanh 


a f 


2 (A y y) 

Ai/ 


- 1 


2(A * - z) , 

a f ( t 1 


)]}• 

)]}■ 


( 3 ) 

( 4 ) 


where af = 2.5, X y = 30 A y /N y , X z = 30 A Z /N Z and A y and A z are respectively the lateral 
and vertical domain sizes. 

The effectiveness of this scheme was tested by applying it to the internal wave system 
created by a point vortex (a multiple- frequency analogue of the vibrating-cylinder, ‘Saint 
Andrew’s cross’, experiments of Mowbray & Rarity [27] (also in [28] and [29]), to verify that 
the internal waves neither re-entered the domain or reflected off the fringe boundaries. The 
energy absorbed by the fringe was calculated during the turbulent wake simulation. The 
strongest activity was detected during the NEQ phase of the simulation, when animations of 
quantities such as dw/dz (cf. figure 12 b) suggest that the waves are generated by individual 
turbulent structures (rather than, say, mean collapse of the wake defect associated with entry 
into the Q2D regime). For the present high-F flow, the internal- wave energy was found to 
be only a small percentage of the overall dissipation (see discussion regarding figure 3b). 

Another computational challenge associated with this flow is defining an appropriate 
domain. The streamwise two-point correlations in figure 2a indicate that the streamwise 
domain size A x is adequate throughout the simulation, and individual structures fit com- 
fortably within the domain (since their streamwise integral scale is much smaller than A x ), 


(«) {b) 




Figure 2. Evolution of (a) streamwise two-point autocorrelations and (6) streamwise spectra 

of streamwise velocity: /o, Nt = 0.2 (in b only); /□, Nt = 3.2; /o, 

Nt — 21; /•, Nt — 136; /♦, Nt = 587. Straight solid line in (6) indicates 

„- 5 / 3 behavior. Spectra En and correlations Rn = Qn/Qn(0) are from averages of 
En and Qu over lateral and vertical planes though the wake centerline in regions for 
which u > 0.5 Udi where u'u' = Eii(k x ) and Qn(r x ) — u'(x)u'(x + r x ). Spectra and 
wavenumbers in ( b ) normalized by U r = 2.09 Udo and L r — h Q . 


even at very late times. The large-scale waviness during the late stages is due to the rel- 
atively low Reynolds number (i.e. narrow range of streamwise scales), and especially the 
periodicity of the large pancake-eddy structures present in the Q2D regime (figure 6). How- 
ever, the streamwise domain size also directly affects the quality of the mean statistics (and 
possibly the dynamics of the merging of the large-scale structures), since it defines the size 
of the ‘large-eddy sample’ (which will decrease as time passes) that enters the streamwise 
mean. Therefore, the choice was made to have a relatively large streamwise extent, with 
A x /h Q = 967 r ~ 302 (cf. Gourlay et aV s [13] A x /h Q = 245, Brucker & Searcher’s [9] 
A x /h 0 ~ 105, and Diamessis et a/.’s [2] A x /h 0 « 54). 

The possibility that the height and especially the width of the wake will outgrow the 
periodic domain is of even greater concern. To avoid this problem, without unduly wasting 
computational resources, the lateral and vertical domain sizes A y and A z are dynamically 
increased at appropriate times during the run, by projecting the periodic solutions onto a 
domain of greater extent in y , and lesser in z, with the same grid spacing. Details of this 
re-projection procedure are given in RCC. The spatial resolution is also reduced at certain 
times during the simulations, as allowed by the monotonically decreasing Reynolds number 
(see table 1). 

The choices for the lateral and vertical domain are tabulated in table 1, with respect to 
both the initial h Q and local mean lateral h y and vertical h z half- widths. The criterion found 
by RCC for the unstratified case, that the lateral and vertical width h must remain less than 
0.17 of the corresponding domain size, is satisfied throughout the simulation. Confirmation 
that stratification does not significantly alter this guideline is provided in figure 3a, where 
the DNS histories are given along with those from a test run identical to the case defined 
in table 1, except that the lateral and vertical domain sizes remain fixed at their values at 
Nt = 28 (corresponding to the end of the NEQ regime). There is good agreement between 
the two runs (including for the maximum Udi width h y and height h z of the mean defect, 
as well as the streamwise u ' , lateral v' and vertical w' velocity fluctuations) until Nt ~ 340, 
at which point h y /A y « 0.15 and h z /A z « 0.12 for the smaller /fixed- domain results. These 
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Nt 

M x 

Ax j T/min 

My 

A7//77 min 

Ay / Hq 

hy/^y 

Ay/hy 

M z 

Az/f] min 

A Z / ho 

hz/ A-z 

A z/h 

0 

3072 

3.8 

256 

1.0 

2tt 

0.159 

0.02 

256 

1.0 

27 r 

0.159 

0.02 

0.87 


7.0 


1.7 


0.155 

0.02 


1.7 


0.157 

0.02 

0.87 



384 


3tt 

0.103 

0.02 

384 


37 r 

0.105 

0.02 

1.12 


6.7 

256 

2.5 


0.112 

0.02 

256 

2.5 


0.110 

0.02 

1.12 



384 


4.57 r 

0.075 

0.02 

384 


4.57T 

0.073 

0.02 

2.86 


2.9 

256 

1.6 


0.131 

0.02 

256 

1.6 


0.120 

0.02 

2.86 



384 


6.75tt 

0.087 

0.02 

384 


6.75tt 

0.080 

0.02 

5.38 


1.8 

256 

1.5 


0.132 

0.02 

256 

1.5 


0.087 

0.03 

5.38 



384 


10.137T 

0.088 

0.02 




0.087 

0.03 

9.72 

2560 

1.4 

256 

1.5 


0.108 

0.02 

192 

1.3 


0.085 

0.04 

9.72 



384 


15.197T 

0.072 

0.02 




0.085 

0.04 

18.6 

2304 

1.0 

256 

1.4 


0.078 

0.03 


0.8 


0.078 

0.05 

18.6 



384 


22.78tt 

0.052 

0.03 




0.078 

0.05 

28.1 

1536 

1.2 

320 

1.4 


0.057 

0.04 

128 

1.0 


0.077 

0.07 

28.1 



480 


34.17tt 

0.038 

0.04 

192 


10.137T 

0.052 

0.07 

35.4 

1024 

1.7 

384 

1.6 


0.040 

0.04 


1.0 


0.051 

0.07 

46.1 

512 

3.3 

256 

2.4 


0.045 

0.06 

128 

1.4 


0.051 

0.10 

66.0 

384 

4.3 

192 

3.0 


0.051 

0.07 

96 

1.8 


0.052 

0.13 

66.0 



256 


45.567T 

0.038 

0.07 

128 


13.57T 

0.039 

0.13 

109 


3.5 

192 

3.4 


0.066 

0.05 

96 

2.0 


0.049 

0.14 

109 



256 


60.75tt 

0.050 

0.05 

128 


18.0t r 

0.037 

0.14 

215 

256 

3.1 

192 

2.6 


0.072 

0.05 

96 

1.5 


0.056 

0.12 

215 



256 


8l7T 

0.054 

0.05 

128 


247T 

0.042 

0.12 

587 


1.5 


1.3 


0.081 

0.03 


0.8 


0.067 

0.08 

841 


1.3 


1.1 


0.094 

0.03 


0.6 


0.079 

0.07 


Table 1. Numerical parameters. The numbers of streamwise M x , lateral M y and vertical M z Fourier-expansion coefficients are two-thirds of the 
corresponding collocation/quadrature points: N x = 3M x /2, N y = 3M y /2 and N z = 3M z /2. Streamwise Ax, lateral Ay and vertical A z grid 
spacings given by A X /N X , A y /N y and A Z /N Z , respectively. Minimum Kolmogorov length scale ? 7 m i n a (z^ 3 / e max ) 1//4 , where e max is the maximum 
of the streamwise-mean rate of dissipation of TKE at each time. The length A x of the streamwise domain is fixed at 967 rh 0 , where h Q is the 
mean wake half- width/height at Nt m 0, while A y and A z increase at times shown (indicated by vertical gaps between rows). Values remain 
unchanged until a new entry is given. 


Nt 




Figure 3. (a) Comparison of results from the present DNS and a fixed-domain test run 
beginning from Nt = 28.1, with A y and A z maintained at 34.177rh 0 and 10.137 rh OJ respec- 
tively, such that the domain is not expanded at Nt = 66, 109 and 215; see table 1): lines, 
present DNS; symbols, fixed-domain test run (o, h y ; •, h z \ ♦, Ud ; o, ( u'u / ) 1 ^ 2 ; ■, (xV) 1 / 2 ; 
□, ((w'w') 1 / 2 , where streamwise (u'u'), lateral (v'v') and vertical (w'w') Reynolds stresses 
are cross- wake averages over the region in the y-z plane where u > 0.05 Ud)- (b) Histories 
from the present DNS of the rate of change of integrated TKE {/c} and net production 

minus dissipation plus buoyancy {V t ke} — {<uke} + {£>tke} (see §3.5 for definitions): , 

d{k}/dt; ♦, {V t ke} - {e t ke}; o, {P t ke} - {etke} + {^tke}- Budget terms in (5) given in units of 
£7 3 L r , where reference velocity U r = 2.09 Udo and length L r — h Q . 


ratios are larger than those produced by the expansions used for the NEQ and Q2D phases 
of the present DNS, which implies the results have not been adversely affected by a too-small 
domain. 

Table 1 also shows the changes in resolution over the course of the simulation. The 
streamwise Ax, lateral Ay and vertical Az grid spacings are small enough multiples of the 
minimum Kolmogrov scale ? 7 m i n throughout the simulation (between 1 and 4, except for a 
brief period at very early times when Ax reaches 7?7 m i n ) that we can safely assume the results 
faithfully capture all the significant spatial scales of the turbulence. The temporal scales 
are resolved by ensuring the Courant-Friedrichs-Lewy (CFL) number does not exceed y/2&. 
These conclusions are supported by the streamwise spectra in figure 2b. The signature of 
the vortex-ring characterize is apparent at Nt = 0.2. The spectra also indicate the moderate 
Reynolds number range considered here, in that ft -5 / 3 behavior is apparent, but only at 
early times, bracketing Nt = 3. A more empirical check of the resolution is provided by 
the histories in figure 36, comparing the sum of the production, dissipation and buoyancy 
terms (the fringe contribution is not included) in the integrated TKE equation (the open 
circles) with its integrated rate of change (solid line) . The good agreement between the two 
points to the accuracy of the present DNS. It also implies that the energy of the internal 
waves generated by the turbulence during the NEQ phase is negligible, since the integral of 
the flux-divergence terms is essentially zero compared to that of the production, dissipation 
and buoyancy. This is in contrast to the low-F case, where departure from the 3D regime 
is accompanied by significant TKE transfer away from the wake (Brucker & Sarkar [9]). 

We note that figure 3 b also reveals the integrated buoyancy term (the difference between 
the solid diamonds and open circles) is negligible except between Nt ^1.5 and 15, and that 
during this time it acts solely as a sink of TKE. This provides an answer to the question 
raised in the introduction, regarding the defining NEQ mechanism, since it indicates that 
the transition from the 3D to the NEQ regime is not accompanied by a net transfer of 
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potential to kinetic energy. 

A final comment regarding the present simulation strategy is warranted: the lateral- and 
vertical-domain expansion is accompanied by a coarsening of the resolution across the mean 
wake profiles u(y, zo) and u(yo, z ), quantified in table 1 by a tendency for A y/h y and Az/h z 
to increase with time. A check was made that the loss of resolution does not compromise the 
shape of the u profiles, by repeating the portion corresponding to the largest A z (> 0.1 h z , 
approximately) with double M z , and confirming that the u results have negligible differences 
compared to those shown below. 


3 Results 

3.1 Flow visualisation 

The turbulence structures associated with the 3D, NEQ and Q2D regimes are shown, respec- 
tively, in figures 4, 5 and 6. The weak stratification used here allows the wake turbulence to 
develop a mature three-dimensional state before buoyancy begins to alter the flow. At the 
time shown in figure 4 (Nt ~ 3) the stratification causes the height to be slightly smaller 
than the width (see also upper right-hand-side inset plot in figure 21 6). However, it is evi- 
dent that there is a substantial range of scales (cf. figure 26), and that the wake is to a good 
approximation still axisymmetric. 

The large-scale lateral oscillation, as well as the flattening and vertical striations, ob- 
served in the NEQ structure (figure 5, for Nt ~ 21) are typical of those found in other 
experiments and simulations [2, 12,30], some of which correspond to lower and higher 
Roo • In contrast to the more isotropic eddies seen at Nt = 3, here the turbulence in figure 5 
contains ‘arrowhead’ structures in the vertical plane through the wake centerline that point 
in the direction of motion of the virtual body associated with the wake (cf. figures 10-12 of 
Diamessis et al. [2]). A visual estimate of the height of the turbulent region (as opposed to 
that of the mean profile) indicates it is similar to, perhaps slightly less than, that for the 3D 
regime (compare figure 46). The change in structure corresponds to enhanced production 
of TKE and profound anisotropy in its rate of dissipation (figure 166). The lateral vorticity 
contours in figure 56 suggest the presence of Kelvin- Helmholtz (K-H) billows at the top 
and bottom of the wake, in the inclined shear layers of the arrowhead structures (e.g. near 
x = 90 h 0 ). Diamessis et al. [2] discuss the manner in which the characteristics of these 
inclined shear layers depend upon Reynolds number, and the associated question of the how 
the layers’ overturning via K-H instabilities can contribute to the generation of turbulence 
in stably stratified environments at high Reynolds number. For the present case, we show 
below that at the time corresponding to figure 5, viscous effects are not insignificant. 

The Q2D regime contains flattened, large-scale, pancake- type motions in horizontal 
planes, such as those shown for Nt ~ 587 in figure 6. There are three separate vortical struc- 
tures that have the two-layered feature that emerged earlier in the simulation (figure 5c), 
which bring to mind Godeford & Cambon’s [11] distinction between two-dimensional and 
two-component turbulence (see also Spedding [12]). These two-layer eddies are similar 
to those observed by Spedding [1] (for Nt > 50 and > 10), by Gourlay et al. [13] 
(Nt > 100, Fqo > 150), by Brucker & Sarkar [9] (Nt = 125, = 8), and by Diamessis et 

al. [2] (Nt = 70, Foq = 4). Figure 6 is also consistent with the vertical layering of lateral 
vorticity fluctuations seen by Spedding [12] (at Nt = 64 for F^ =4). In light of the vortex- 
ring initialization used here, which injects structures very different from the helical modes 
associated with the towed-sphere wake (cf. figure 5a of RCC), the present results add to the 
growing body of evidence that the Q2D eddies are unrelated to the type of instability by 
which the wake became turbulent (see also [2, 13]). 
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Figure 4. Instantaneous flow at Nt = 3.2 (Nt = Nt — Nt vo = 2.5, end of 3D regime; 
virtual origin t vo defined in figure 7). (a) Top and ( b ) side views of isosurfaces of vorticity 
magnitude (|u;| = 0.38 U r /L r ). Mean-flow direction defined as if a virtual wake-creating 
body has moved from left to right through the sub-domain shown. Views shown are 8% of 
streamwise, and 47% of lateral and vertical, domain. 
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Figure 5. Instantaneous flow at Nt = 20.9 (Nt = Nt — Nt vo — 20.2, NEQ regime; virtual 
origin t vo defined in figure 7). (a) Top and (c) side views of isosurfaces of vorticity magnitude 
(|u;| = 9.7 x 10 ~ 4 U r /L r ) . ( b ) Contours of lateral vorticity uo y in vertical plane through 
centerline, in units of 10 2 U r /L r . Motion of virtual wake-creating body is from left to right. 
Views shown represent subregion of full domain (table 1). 
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Figure 6. Instantaneous flow at Nt = 586.9 (Nt = Nt — Nt vo — 586.2, Q2D regime; 
virtual origin t vo defined in figure 7). (a) Contours of vertical vorticity uj z in horizontal 
plane through centerline, in units of 10 A U r /L r ). ( b ) Top and (c) side views of isosurfaces of 
vorticity magnitude (|ca| = 3.1 x 10 -7 U r /L r ). Motion of virtual wake-creating body is from 
left to right. Views shown represent full streamwise, but subregions of lateral and vertical, 
domain (table 1). 
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3.2 Mean defect history 

The history of the maximum mean velocity defect Ud is shown in figure la and b. The Ud 
variation is smoother than that typically found from simulations, due to the reduced sta- 
tistical oscillations associated with the larger-than-usual streamwise domain used here (the 
same can be said for the mean- wake height and width histories in figure 7c, d). In contrast 
to the unstratified- wake DNS of RCC, for which the difference between the maximum and 
the centroid of the mean defect increased with time, and becomes significant at later stages, 
here there is essentially no difference between the two, for the time period considered. There 
is however a noticeable difference in the location in the y-z plane at which the two occur; 
as in RCC, we use the location of the centroid of the mean-velocity distribution u(y,z) to 
define the current /local wake centerline coordinates ( y c ,z c ). 

Figure 7b also shows histories of the streamwise u' , lateral v' and vertical w' velocity 
fluctuations at the wake centerline (symbols). These, along with the Ud history, reveal that 
it is not until after Nt ~ 1 that the wake turbulence is fully developed. We therefore view 
the early phase, from Nt = 0 until Nt = Nt vo 1, as effectively a precursor simulation. The 

3 / 2 

exact virtual-time origin t vo is defined by the unstratified, 3D similarity diagnostic U d 7 
shown in figure 7a, from which we deduce Nt vo = 0.645 (chain-double-dotted line). In order 
to remove some of the uncertainty associated with comparing the present results during the 
3D and NEQ regimes with those from experiments and other simulations, the histories that 
follow will be presented as functions of Nt = Nt—Nt vo = Nt— 0.645. Figure 7 a also indicates 
that Ud continues to exhibit pure-3D behavior until Nt = Nt — Nt vo « 4.8 — 0.645 « 4.2, 
since this is when the similarity diagnostic first deviates from the linear variation. 

The early wake is well represented by the canonical Ud ~ t -2 / 3 history despite not being 
fully self similar during this time: although fully turbulent and essentially unstratified and 
axisymmetric (figure 8a), the wake’s non-dimensional Reynolds-shear-stress profile in cylin- 
drical coordinates —v' x v' r /U d (not shown) is not solely a function of the non-dimensional 
radius r/h (i.e. with no time dependence). This is consistent with the evolution of the un- 
stratified, vortex-ring-initialized, Case VR wake described in RCC, for which (non-universal) 
self-similarity did not appear until after the time corresponding to Nt « 11 for the present 
simulation - and universal self- similarity not until Nt > 70. (For the present DNS, the 
relationship between Nt and the non-dimensional time used in RCC is t/t* = 30Nt.) In 
fact, in the range corresponding to 1 < Nt <11, the Ud history for Case VR decays as £ _1 , 
because of the unique turbulence distribution induced by the Case VR rings (which produce 
an eddy viscosity profile that is temporally constant and radially uniform; see figures 8 a and 
9 b of RCC). Here the ring initialization, which used somewhat larger streamwise spacing, 
2.1 h Q rather than Case VR’s 1.85 h Q , is such that Ud ~ t~ 2 ^ 3 is a better representation than 
t~ l for the 3D regime. These observations illustrate the profound impact the initial wake 
structure can have upon its development. 

As observed elsewhere (e.g. in Spedding’s [1] experiments and Gourlay et a/.’s [13] 

simulation), the magnitude of the rate of decay \dUd/dt\ slows markedly for Nt > 5, as 
the flow enters the NEQ phase, such that Ud ~ t~ n with n < 1/4. In this region, the Ud 
variation is well represented by the interpolant 0.25 (AT + l) -1 / 4 . (It should be mentioned 
that power laws with slightly different exponents could also reproduce the Ud decay to good 
accuracy, given the flexibility afforded by the use of a finite virtual origin for time.) The 
time at which the 3D and NEQ power laws intersect, Nt = Nt — Nt vo « 6 — 0.645 « 5.5, 
can be viewed as an estimate of the boundary between the two regimes for the present flow. 
The generality of this finding will be discussed in §3.6. 

The Ud history for the variable- viscosity run (Case R1500; see § 2) begins to significantly 
differ from the constant-viscosity/variable-i^ result after Nt « 100 (compare dotted and 
solid lines in inset plot in figure 7b), which implies non-inertial effects are important after 
this time. The NEQ period revealed by the evolution of Ud begins at Nt ~ .5.5 and 
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Figure 7. Histories of (a, b) maximum mean velocity Ud , (c) mean height l z and vertical 

spreading rate /3 Z and (d) mean width £ y and lateral spreading rate /3 y : , present 

DNS (with £ z = r* and £ y = r* in c and d, respectively); □ (in b only), 0.1 (u'u'^J 2 /U a 
(arbitrary factor of 0.1 is to aid visibility); A (in b only), 0.1 /Ud\ ▼ (in b only), 

0.1 (w'w')l^ 2 /Ud] (in c and d), h z and h y \ (in c and d), r' z and r’ y \ o (in d only), 

L' a /r y (cf. Meunier & Spedding [17]); • (in d only), L"/r* (cf. Meunier & Spedding [17]); ♦, 
(in c) p z = (l/U d )dr* z /dt and (in d) P y = (l/U d )dr*/dt; o, (in c) /3 Z = (l/U d )dr* z /dt and (in 

d) f3 y = (l/Ud)dr*/dt for Case R1500 (variable- viscocity (C/^r */ v > 1500) run); , 

curve fits for (a, b) Ud , (c) r* and (d) r*, where in (a) ( Ud/Ud(t))~ 3 / 2 = 3.1(At— 0.645), in ( b ) 
U d /U do = 0.47 (At — 0.645) _2 / 3 for At < 10, U d /U do = 0.25(At+ 1.0)- 1 / 4 for 3 < Nt < 100, 
U d /U do = 2.05(At-10)-°’ 85 for At > 40, while in (c) and (d) r*/d 0 = rj/d 0 = 1.23( At) 1 / 3 
for At < 10 where At = At — Nt vo and Nt vo = 0.645 (see figure 215). © (in b only), time at 

which Case R1500 first reaches its minimum Reynolds number, U d r* / v — 1500; (inset 

plots in 5, c, d only), Case R1500 (with £ z m r* and i y — r* in c and d). At At = 0, 
r t ~ r y — r o ~ 0.865d o . The U d (t) interpolants intersect at At « 6 and 45. 
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follows (or remains just above) the NEQ power-law curve fit until about Nt = 30 and 
Nt = 50 respectively for the constant- and variable- viscosity (i.e. lower and higher Reynolds 
number) cases. The effect of viscosity on the NEQ duration is even more apparent in 
the height histories (figure 7c). The Reynolds numbers for both cases are thus too low to 
allow an extended NEQ region - and that the period during which significant turbulence 
is ‘regenerated’ by Kelvin-Helmholtz instabilities (Diamessis et al. 2011), if it exists at all, 
will be quite brief (cf. figure 5 b). (see instantaneous lateral vorticity contours at Nt = 20 in 
figure 5 of Redford et al. 2014). The Reynolds- number dependence of the duration of the 
NEQ regime will be considered further in §3.6. For now we note that the Ud history seen 
here during the final stage of the NEQ phase is more consistent with the £ -1 / 2 behavior of the 
‘buoyancy-controlled’ subregion between NEQ and Q2D proposed by Meunier et al. (2006) 
than the Ud plateau of the ‘NEQ 2 ’ subregion identified by Brucker & Sarkar (2010). 

A £ -3 / 4 power law for Ud could be assumed for the Q2D regime (cf. [3,25]), which would 
agree well with the present data for Nt > 80 if a virtual origin of Nt = 30 were used. 
(Recall the earlier comment about the interplay between the power law’s exponent and 
virtual origin.) However, a better fit is provided by £ -0-85 for Nt > 60, with a smaller virtual 
origin (Nt — 10) (chain-double-dotted curve in figure 7a). The —0.85 exponent is also more 
consistent with conservation of mean momentum, when coupled with the behavior of the 
vertical and lateral defect dimensions (see below). The NEQ/Q2D boundary defined by the 
intersection of the corresponding power-law interpolants is Nt ~ 55; the t -0-85 interpolant 
closely corresponds to the Ud variation until the end of the (primary) simulation, Nt ~ 840. 
Note that the final time is not large enough for a purely viscous 3D regime to appear (such 
that Ud ~ £ -1 ; see Meunier et al. [3]). In fact, the TKE remains a significant fraction of 
C/J for all the times considered here (figure 15a inset plot). Nevertheless, the importance of 
viscous effects during the Q2D phase is underscored by the significant deviation for Nt > 90 
of the Ud histories from the constant- and variable- viscosity flows (figure 7b inset plot). 

The height histories are shown in figure 7c. Both the ‘half-radius’ height h z (i.e. half the 
distance between the locations at which the mean velocity in the vertical plane through the 
centerline, u(z, y c ), is 0.5 of Ud ; not a best-fit to an assumed Gaussian profile), and an inte- 
gral measure r* are shown; the latter is defined as r* = [(2 In 2 /Ud) / 0 °° u(r z , 2 / c )'r^d.r*^] 1//2 , 
where r 2 = (z — z c ) 2 and (y c ,z c ) are the coordinates of the local wake centerline. The 
integral height has been normalized such that it is equal to h z when u(z, y c ) is Gaussian. 
Consequently, the greater the difference between h z (chain-dotted line) and rj (solid) in 
figure 7c, the more non-Gaussian the vertical mean velocity profile (figure 9a). 

During the 3D phase, when the wake is effectively unstratified, its height and width 
grow in proportion to t 1 / 3 , and its non-dimensional mean growth/spreading rate /3 Z = 
(l/Ud)dr*/dt is approximately constant. This can be seen in figure 7c, for Nt < 3 (chain- 
dotted curve and symbols, which respectively trace a t 1 / 3 interpolant for rj and the /3 Z 
variation) . Perhaps unsurprisingly, given the flattening and widening effect the stratification 
has on the mean defect, the height /width interpolant deviates from the DNS sooner (near 
Nt — Nt vo = 2) than the corresponding Ud power-law interpolant does (Nt ~ 4.5). 

Despite the wake not becoming fully self-similar, the spreading rate during the unstrat- 
ified regime remains fixed near 0.03, which is roughly half that observed at much later 
times (corresponding to Nt > 70) for the self-similar growth of RCC’s unstratified Case VR 
wake. The streamwise velocity fluctuations at the centerline are also weaker here during 
the 3D regime, than for Case VR during its (first, non-uni versal) fully self-similar phase, 
with (u'u' )l /2 /Ud ~ 0.2 here (open squares in figure 7b) compared to approximately 0.3 for 
Case VR. 

The difference between the integral height for the constant- and variable- viscosity runs 
(solid and dotted curves, respectively, shown in inset plot), and that of the non-dimensional 
mean growth/spreading rate /3 Z (solid and open diamonds) - namely that the region of zero 
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slope (and thus zero (3 Z ) lasts longer for larger Rd ( Nt 70 rather than 30) - corroborate 
Meunier et aids [3] prediction, and Diamessis et aids [2] observation, that the duration of 
the NEQ regime increases with increasing Reynolds number. 

Figure 7c also shows the differences in the evolution of the height of the mean defect 
and vertical extent of the wake turbulence, which become apparent when one compares 
histories of h z (chain-dotted curve) or r* (solid) with that of r z (dashed), a TKE based 

integral length, defined as r' z = ({ k (z — z c ) 2 }/{k}) where k = \u\u\. The mean defect 
does indeed ‘collapse’ early in the NEQ regime, in the sense that both h z and r* decrease 
with time from Nt ~ 8 until about Nt = 15 - although the reduction in the integral height 
is much weaker than that exhibited by the u = 0.5/7^ measure h z . As found by Gourlay et 
al. [13], the turbulent region continues to spread vertically during the NEQ regime, when 
dr*/d£ ~ 0, before an even more dramatic collapse occurs, with dr z /dt < 0 (much larger 
than the negative dh z /dt and dr*/d£ seen earlier) between Nt ~ 20 and 100. This behavior 
was also observed for the self-propelled case by Lin & Pao [10] in their flow- visualisation 
experiments, and in Brucker &; Sarkar’s [9] DNS of self-propelled and towed wakes. 

After the wake enters the Q2D regime, r' z begins to grow again, as do h z and r*. All 
three heights then converge, as the vertical spreading becomes a purely viscous phenomenon, 
such that h z « r* « r' z [13] and each grow as t 1 / 2 [3]. This can be seen in figure 8, which 
includes contours of the mean defect (lines) and TKE (colour) at times representing the 
three regimes. The tendency for the height of the turbulent region to occupy a smaller and 
smaller fraction of the mean defect as time passes is clear. (The changes to the TKE budget 
corresponding to the change in shape of the TKE profiles - for example, the growth of the 
‘flattened, double-pronged’ structure in the y-z plane that begins to appear during the NEQ 
phase - will be considered in § 3.5.) The non-dimensional vertical growth rate /3 Z is therefore 
positive but much smaller than during the 3D regime, when the spreading is effected by the 
turbulence. 

Due to the conservation of mean streamwise momentum, the product of the mean height, 
width and maximum velocity will remain approximately constant through the wake’s life- 
time. (It would be exactly constant for a pure 2D-Gaussian mean defect. The Udh z h y 
product for the primary simulation remains within ±16% of its value at Nt = 1 for the 
full 1 < Nt < 840 range.) The variations of the power-law exponents for Ud and h z during 
the 3D, NEQ and Q2D regimes are thus expected to lead to variations for the mean- width 
history (which would conflict with Diamessis et a/.’s [2] claim that h y ~ t 1 / 3 for all three 

regimes). For example, the combination of monotonically decreasing Ud and dr^/dt — »• 0 
requires the mean width to grow faster than it does during the 3D phase as the wake en- 
ters the NEQ regime. Figure 7 d illustrates the enhanced lateral growth, compared to the 
3D power law (chain-double-dotted curve), in terms of both the mean half-width/radius 
h y (chained-dotted curve) and the integral width r* & [(2 In2/Z7d) u(z c ,r y )r y dr y ] 1 / 2 , 
with Ty = {y — y c ) 2 (solid curve), for 1.5 < Nt <4.5. Thereafter, the rate of spreading 
slows, falling from a maximum near Nt = 5 to a local near-zero minimum near Nt = 20 
(see open and closed diamonds in figure 7 d, which show the history of the non-dimensional 
lateral spreading rate /3 y = (l/Ud)dr*/dt). Beginning from Nt ~ 25, as the wake leaves 
the NEQ and enters the Q2D regime (and Ud(t) begins to increase its rate of decay), the 
lateral growth increases rapidly until it reaches a value comparable to that observed during 
the unstratified 3D period, with r* ~ t 0-35 . This is similar to the h y ~ t 1 / 3 proposed by 
Diamessis et al. [2], but differs from the t 1 / 4 included in Meunier et aids [3] model, and the 
£°- 23 observed in Brucker & Sarkar’s [9] DNS, for the Q2D range. The difference between 
the n = 0.35 and 0.25 exponents is shown by the thin-solid lines in figure 7 d, which supports 
our preference for the larger value - which in turn explains (in light of the viscous n = 1/2 
growth of h z , and conservation of momentum) our use of n = —0.85 to fit the Ud(t) history 
(figure 7b). (Meunier et aids [3] exponent combination of (—3/4, 1/2, 1/4) for ( Ud,h z ,h y ) 
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(a) Nt — 1.2 ( b ) Nt = 2.5 



(d) M = 82 



(e) M = 586 
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Figure 8. Mean-defect and TKE contours at: (a) Nt = 1.2, ( b ) Nt = 2.5, (c) Nt = 18, 
(d) Nt = 82 and (e) Nt = 586. TKE (shaded contours) normalized by local C/J. Solid-line 
contours correspond to u — 0.05 Ud and 0.5 Ud- 
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also satisfies the constant-momentum constraint.) Despite the profound qualitative differ- 
ences in the turbulence structure in the 3D and Q2D regimes, the rate at which the mean 
defect spreads laterally is about the same. 

The width of the turbulence, like that of the mean defect, grows continuously throughout 
the lifetime of the wake. This can be seen by the history of the TKE-based integral width 
r ' y — ({k (v ~ yc) 2 }/{k}) in figure 7 d (dashed line). In contrast to the height of the 
turbulent region, its width remains significantly larger than the width of the mean defect, 
with r' y /hy > 1 for each of the 3D, NEQ and Q2D periods (cf. Gourlay et al. [13]). (We 
cannot say whether or not the waviness in r' y (t) for Nt > 20 is physical.) The closed and open 
circles in figure 7 d respectively correspond to the ratios L'^/r* and L"/r*, where L' a and L" a 
are, respectively, the distance from the wake centerline to the maximum values of the lateral 
Reynolds shear stress —u'v' and the horizontal-plane velocity fluctuations (u'u' + v'v') 1 / 2 . 
Integral version of these scales are used, based on the profiles presented in Meunier & 
Spedding [17], with L' a = \u'v'\dy/u 0 2 and L" a = /^(uV + v'v'f^dy / (2 'k) 1 ^ 2 Uq, 

where Uq = exp(0.5)(?/V) max , — 0.5 exp (0.5) (id id + i/i/)max, and (u f v f ) max and (u'u' + 

ijVjmax are the maxima over y of \u'v'\ and (u'u' + v'v') 1 / 2 , respectively. 

In contrast to r^/r*, which is always significantly greater than one, the shear-stress width 
ratio L' a /r* remain close to unity for the entire run, slightly less before Nt « 10, slightly 
more thereafter, while for the horizontal TKE measure 0.7 < L"/r* < 1 throughout. The 
L' a ~ r* and L " ~ r* behavior, which is similar to that exhibited by the experimental data 
for various bluff bodies given in Meunier &; Spedding [17], implies a type of self-similarity 
for the lateral profiles, especially for — u'v', that survives as the importance of stratification 
and viscous effects grows, and the flow moves through the 3D, NEQ and Q2D regimes. 

The decoupling of the vertical and lateral spreading processes, and the associated greater 
importance for the latter of the mixing due to turbulence relative to molecular viscosity, 
is illustrated by the much closer agreement of the mean width r* from the constant- and 
variable- viscosity flows (figure 7 d inset) compared to the corresponding r* histories (figure 7c 
inset). 

Figure 9 contains vertical (9a) and lateral (96) profiles of streamwise-mean velocity 
u, normalized by the local maximum Ud and height h z or width h y of the defect. The 
better collapse and smoother profiles seen here, compared to those from Gourlay et a/.’s 
[13] DNS and Diamessis et a/.’s [2] LES, are thought to be a symptom of the better 

averaging sample afforded by the present larger streamwise domain. Included in these plots 
are results from times corresponding to the structures shown in figures 4 (dashed curve; 3D 
regime), 5 (chain-dashed; NEQ) and 6 (long-dashed; Q2D). (The dotted profile is the mean 
induced by the vortex ring initialization.) Given the qualitative differences in the 3D, NEQ 
and Q2D structures, the collapse at these times is striking, as is the agreement with the 
Gaussian idealization (open symbols), particularly near the centerline. At the earlier times, 
the non-Gaussian behavior away from the centerline may be an artifact of the vortex-ring 
initialization, and/or due to attenuation of the eddy viscosity, associated with corrugations 
of the interface between the vortical and irrotational regions [31] (figure 4). 

For the lateral profiles, which tend to be more consistently Gaussian than the vertical 
ones, the near-centerline agreement is similar to that found by Spedding et al. [25]. (The 
unphysical lateral asymmetry at Nt = 137 is evidence of finite eddy sample at this time. 
Curiously, the u(z c ,y) profile at Nt = 586, which is the result of averaging over an even 
smaller number of eddy structures than are present at Nt = 137, is nearly exactly symmetric; 
see long-dashed curve in figure 96.) The vertical profiles, on the other hand, can at times 
exhibit significant non-Gaussian features, across the entire wake. One might suspect the 
‘extra peakiness’ near the centerline during the early stages of the Q2D regime (Nt = 82 
and 137) is due to the coarsening of the vertical resolution across the mean defect (i.e. 
A z/h z = 0.14; cf. table 1) mentioned at the end of §2. However, this possibility was ruled 
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(a) (b) 




( z z c )/h z (y y^)fhy 

Figure 9. Profile of mean streamwise velocity in (a) vertical and ( b ) lateral planes through 

the wake centerline: , Nt = 0; , Nt — 1.2 , Nt = 2.5; , Nt — 

20; , Nt = 82; , Nt = 137 , Nt = 586; o, Gaussian idealisation 

u/Ud = exp(— ln(2) s 2 ), where s is either ( y — y c )/h y or (z — z c )/h z . Coordinates (yo,zo) 
and (y c , z c ) respectively define the initial and current location of the wake centerline, given 
by the centroid of the mean velocity defect. 


out by repeating this portion of the simulation with increased M z , such that A z/h z = 0.07, 
resulting in u(z, y c ) variations that are nearly indistinguishable to those shown in figure 9a. 

The mechanisms responsible for the evolving shape of the mean defect revealed in figure 9, 
as well as its magnitude, height and width (figure 7), can be further understood by examining 
the corresponding evolution of the terms in the mean- momentum transport equation. But 
before attending to that task, we consider the history and budget of mean kinetic energy. 

3.3 Integrated mean kinetic energy history and budget 

The influence of the stably stratified background on the three components of mean kinetic 
energy (MKE) is shown in figure 10a. Histories of the (cross-domain) area-integrated quan- 
tities {K aa } = {^u a u a } for a = 1,2,3 (no sum on a) are presented. Only during the 
NEQ regime are the lateral and vertical energies non-negligible compared to the streamwise 
component, and even then remain at least 1.5 orders of magnitude smaller than {K\ i} at all 
times. Nevertheless, all three components reveal discernible buoyancy effects. These can be 
deduced from figure 105, which presents histories of the dominant terms in the integrated 
MKE transport equation: 


d {Kgg} 
d t 


where 


— {Paajj }’ {Dcxatjj } T {F aQ }, OL — (1,2,3), 
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is the mean-to-turbulence transfer term (i.e. the negative rate of TKE production), 

n _ du a du a 
Uaajj - V d d 

is the mean dissipation rate, and 

Baa — QP^a Poo ^mke = 9P ^ / Poo 
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Figure 10. (a) Components of integrated mean kinetic energy {K aa } = {^u a u a } (no sum 
on a). Values normalized by MKE at Nt — 0, {Ka} 0 m {Ku} 0 + {^ 22)0 + {^ 33 }o- (6) 
Dominant terms in budget of MKE components, in units of V 3 L r , where U r = 2.09/7^ o and 
L r = h Q . (These arbitrary scales are implicitly defined by the non-dimensionalization used 
to prescribe the radius, spacing and circulation of the vortex rings in the initial-condition 
field.) Note signs of {F 2222 } and {B m ke} = {^ 33 } (open symbols), which (unlike all the 
other terms) both act as sources in their respective +d{K aa }/dt balances. 


is the mean buoyancy transfer (which appears with opposite sign in the transport equa- 
tion for the mean potential energy \gp~p/ Poo |T|). Note that the integral of both the ad- 
vection A aa jj = —UjdK aa /dxj and flux-divergence/transfer T aa jj = —du a u' a u'-/dxj + 
v d[u a du a /dxj\/dxj terms are zero, {A aa jj} = {T aa jj} = 0. For the present flow, for 
which d( )/dxi = d( )/dx = 0, -P aa jj = -P aa 22 - Paa 33 and D aajj = D aa 22 + D aa 33 - 

Regarding the K\\ balance, the end of the 3D period is first evident in the dissipation 
terms {^ 1122 } and {D 1133 }, in that the mean vertical gradients enhanced by the vertical 
flattening of the wake overwhelm the lateral variations, thereby arresting the reduction of 
{D 1133 } (compare single- and double-dot-dashed curves in figure 10b). This begins near 
Nt = 1.5; however, the dissipation terms are too small, relative to the TKE production, to 
have a controlling influence on the {K\i} history. Shortly thereafter, at Nt ~ 2.5 (when 
the Ud history first departs from the t -2 / 3 interplant for the 3D regime; see figure 76), the 
curvature in the {K n} history changes, as a result of near-identical reductions of {P 1122 } and 
{.Pi 133 }, the MKE-to-TKE transfer terms respectively associated with the lateral —u'v' and 
vertical —u'w' shear stresses. Buoyancy has begun to affect the structure of the turbulence, 
such that —u'w' becomes a smaller and smaller fraction of the TKE for Nt > 2.5, while 
—u'v'fk remains nearly constant until well into the NEQ regime (see figure 156). (We shall 
see below (figures 15a and 16) that the TKE begins to increase its rate of decay \{dk/dt}\ 
near Nt = 2.5, as the result of weakening TKE production and growing buoyancy transfer 
(acting as a TKE sink), relative to the local rate of TKE dissipation.) Combined with 
buoyancy- induced alterations of the mean flow - which reduce (compared to the unstratified 
counterpart) the lateral shear du/dy and enhance the vertical shear du/dz (see §3.4 and 
figure 11 ) - this leads to the comparable reductions of — P 1122 = u'v' du/dy and — -P 1133 = 
u'w' du/dz seen in figure 106 (see also figure 166). The result is a less-rapid decay of {K\ 1 } 
(and presumably Ud), as the MKE-to-TKE transfer is diminished, as the wake enters its 
NEQ phase. 

At Nt ~ 5.5, the time at which the 3D and NEQ interpolants for Ud intersect (figure 76), 
figure 106 implies that the processes begun at Nt ^2.5 carry on unabated, until Nt ~ 9 
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(somewhat after the time, Nt ~ 6 . 5 , at which the Ud history first begins to track with the 
NEQ interpolant; figure 75 ), when the buoyancy drives the vertical velocity fluctuations, and 
thus P1133, to zero, leaving P1122 (after Nt ~ 15 ) to account for all the the MKE-to-TKE 
transfer from the streamwise component. (Recall that for this flow the nominal duration of 
the NEQ regime, based on the times at which the 3 D, NEQ and Q 2 D power-law interpolants 
for Ud intersect, is 5.5 < Nt < 45 , while its ‘formal’ duration, based on the period when Ud 
tracks closely with the t -0 25 interpolant, is 6.5 < Nt < 30 .) Dommermuth et al. (2002) and 
Brucker & Sarkar ( 2010 ) have highlighted the correlation between vanishing Pi 133 and the 
occurrence of the NEQ regime; the present results indicate that its inception is triggered 
by a joint reduction of P1122 and P1133, which begins near Nt = 2 . 5 . We note that the 
Pi 122 and Pi 133 histories do not deviate from each other until after the magnitude Ud and 
especially the height h z of the mean defect begin to display characteristics associated with 
NEQ behavior (figures 75 and c), and that this (Nt ~ 9 ) deviation is not only the result 
of the rapid reduction of P1133, but also a consequence of P1122 beginning to ‘wake up’ at 
about the same time (which mitigates the \dUd/dt\ reduction). The attenuation of the rate 
of decay (and eventual change of sign, between Nt « 25 and 60 ) of d{Pn22}/dt can be 
attributed to another buoyancy- forced alteration of the turbulence structure, such that the 
magnitude of the lateral shear stress (figure 155 ) becomes large enough to overcome the 
reduction in the mean lateral shear du/dy and thus produces enhanced Pi 122 • 

By the time Ud departs from the t -1 / 4 NEQ interpolant, near Nt = 30 (figure 75 ), 
the rate of decay of the streamwise MKE, d{Ku}/dt, has returned to levels comparable 
to its pre-NEQ value. This increase is initially due to comparable contributions from the 
lateral MKE-to-TKE transfer P1122 and vertical-gradient dissipation P1133. The ‘bump’ 
observed in the {K\ 1} history between Nt ~ 25 and 100 (see also the Ud history over 
this range, in figure 75 ) corresponds to the period of Pi 122 enhancement, during which 
{P1122}/ {P1133} > 1 , bridging the NEQ-to-Q 2 D transition. (We shall see in § 3.5 that this 
period includes the time, 25 < Nt < 90 , over which the net rate TKE production exceeds 
that of the net TKE dissipation; see figure 16 .) The slowing of the decay of P1122 appears to 
be an integral part of the NEQ regime, with its amplification marking the time at which the 
wake begins to move towards its Q2D phase. For Nt > 150 or so, the differences between 
the lateral TKE production (which has begun to fall again) and vertical-gradient dissipation 
are statistically insignificant. 

The {K\ 1} history thus reflects alterations to both the mean flow (increased du/dz , am- 
plifying the mean dissipation and modulating the MKE-to-TKE transfer; decreased du/dy , 
also modulating the MKE-to-TKE transfer) and turbulence (reduced —u'w' and enhanced 
—u'v', respectively throttling and augmenting the MKE-to-TKE transfer terms Pi 133 and 
Pi 122)- Two further effects of buoyancy can be observed in the vertical and lateral compo- 
nents of the MKE budget. The first is direct transfer via net positive B m ke = — g~pw/ poo, 
which supplies energy from the MPE to the vertical component (open diamonds in figure 105 , 
which for illustration purposes trace positive {B m k e } on the —d{Kss}/dt plot). Note how- 
ever, that in terms of its effect on the mean defect (since k e it is not a direct source 
of Kn), this can be viewed solely as a symptom of the mean-flow collapse, dw/dz < 0, 
associated with the indirect/modulating influences cited above. 

The second further-buoyancy effect visible in the non-RTn budgets in figure 10 is manifest 
by another alteration of the TKE production terms, again due to negative dw/dz < 0 and 
(due to continuity) lateral expansion (dv/dy > 0) of the mean flow in the wake core, such 
that — P3333 = w'w' dw/dz < 0 and — P2222 = v'v'dv/dv > 0 . The latter ‘anti-TKE- 
production’ is responsible for the growth of {PT22 } , the former mitigating the increase in 
{K33} caused by {B m k e }, such that {K33} remains much smaller than {^22}- (Although 
P2222 is a net sink of TKE during the NEQ regime, it is not an important one, as can be 
seen in figure 165 of § 3 . 5 ; P3333, is also insignificant, except within the 40 < Nt < 100 range, 
when it acts as a TKE source, of magnitude of a few percent of the TKE dissipation rate.) 
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Figure 11. Terms in the mean streamwise momentum budget at Nt = 2.5: (a) rate of change 
du/dt ; ( b ) advection —vdu/dy — wdu/dz ; (c) Reynolds-stress divergence —du'v'/dy — 
du'w'/dz ; (d) viscous-stress divergence v(d 2 u/dy 2 + d 2 u/dz 2 ). Terms given in units of 
U 2 /L r , where U r = 2.09 Udo and L r = h Q . 


Having catalogued the above buoyancy effects, the question now arises of how they 
combine to modify the maximum velocity and y-z spatial structure of the mean defect (i.e. 
lid, h z and h y ) during the 3D, NEQ and Q2D regimes. To address this question, we turn 
our attention to the budget of mean streamwise momentum. 

3.4 RANS budget 

In this section we examine the terms in the budget of mean streamwise momentum, and 
how stratification modifies their spatial structure in vertical (y-z) planes normal to the wake 
axis, to better understand the mechanisms responsible for the decay and spread of the mean 
defect during the 3D, NEQ and Q2D regimes. For the present time-dependent parallel flow 
the streamwise Reynolds-averaged Navier-Stokes (RANS) equation can be written 

du ( _du _du\ f du'v' du'w' 

dt \ V dy W dz / + \ dy dz 

advection Reynolds-stress divergence viscous-stress divergence 

During the 3D phase, the reduction in Ud and increase in the wake height and width 
- i.e. the negative du/dt near the centerline and the positive du/dt at the wake edge - 


+ 


d 2 u d 2 u 
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dy : 


dz 2 
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is driven by the divergence of the Reynolds stress, —du'v'/dy — du'w' jdz. This can be 
inferred from figure 11, which shows the cross- wake plane distribution of du/dt and how 
the advection, turbulence and viscous terms contribute to it at Nt = 2.5. At this time, the 
axisymmetric/3D nature of the flow has been slightly altered by the buoyancy, as the flow 
moves towards the NEQ regime. Although the relative strengths of the vertical and lateral 
components of the Reynolds-stress divergence remain comparable at the edges of the wake, 
the flattening and widening of the region of wake turbulence caused by the buoyant forcing 
is apparent (figure 11c). The reduced rate of vertical spreading is due to vertical advection, 
with —wdu/dz < 0 above and below the centerline (dark regions in figure 115), while the 
lateral spreading is enhanced by the corresponding advection term, in that —vdu/dy > 0 
at the left and right of the defect (light regions in figure 115). 

Within the NEQ regime the story is more complex. Figure 12a shows that at Nt = 
18, when dh z /dt ~ dr*/d£ « 0, the regions at the top and bottom of the wake contain 
negligible du/dt , while u continues to grow at the sides. Significant internal wave activity, 
in the k x — 0 (i.e. streamwise mean) mode, is generated by the wake, revealed by the 
vertical mean momentum flux divergence, —duw/dz (figure 125). However, the waves do 
not make a meaningful contribution to the net momentum (or the kinetic energy), since 
—duw/dz — duv/dy = —wdu/dz — vdu/dy is nearly identically zero outside the wake 
(figures 12c and d). The (slower, compared to 3D) rate of decrease of Ud(t) is now driven 
by lateral Reynolds-stress divergence —du'v'/dy and vertical advection — wdu/dz , while 
the lateral component —vdu/dz retards the lateral spreading near the centerline plane but 
increases it further above and below (figures 12 d). The Reynolds stresses, especially the 
vertical component — u'w ', are profoundly affected by the buoyancy (which drives w' — > 0; 
see figure 155), such that the importance of their divergence has been reduced to levels 
comparable to that of the vertical viscous stresses. (Contrast figures 12e, / and g). That 
viscous effects are important within the NEQ regime was illustrated by the difference in the 
defect histories for the primary simulation and Case R1500 (figure 7). The lateral advection 
and the vertical viscous terms (which tend to mirror and counteract each other) reveal a 
‘flattened-X’ pattern in the mean velocity that is not apparent in the vertical u(z, y c ) and 
lateral u(z c , y) profiles in figure 9. The lateral viscous-stress divergence (not shown) remains 
negligible. 

The difference between the duw/dz contours in figures 125 (Nt = 18) and 135 (Nt — 82) 
displays the tendency for the strength of the internal waves to fade as the flow passes from 
the NEQ to Q2D regime. Once the wake enters the Q2D phase, the core deceleration is 
dominated by the lateral Reynolds-stress divergence, and to a lesser extent by the verti- 
cal advection and vertical viscous-stress terms. This can be seen in figures 13e, c and /, 
respectively. The vertical Reynolds-stress —du'w'/dz and lateral viscous-stress i sd 2 u/dy 2 
terms are now negligible. We also note that the multi-layered structure first seen in the 
instantaneous vorticity contours earlier in the NEQ range, at Nt « 20 (figure 5), has now 
matured to the point that it affects the RANS budget, leaving its signature in the terms 
involving spatial gradients of the mean velocity, such as —vdu/dy (figure 13d). At the top 
and bottom edges of the wake, the vertical viscous-stress divergence is now solely responsible 
for the small net positive du/dt here - a fact foreshadowed by the t 1 / 2 growth of the mean 
height observed in figure 7c. In contrast, the (faster) lateral spreading is due to the lateral 
advection (figure 13d) and especially the lateral Reynolds-stress term (figure 13e). 

As time passes, the importance of the —du'v'/dy and vd 2 u/dz 2 terms continues to 
grow relative to the other ones, until they are the only significant contributors to the RANS 
budget (cf. Meunier et al. [3]). At Nt = 586 (figure 14), both decelerate the core region by 
about the same amount, but have very different effects on the spreading of the defect: the 
mean acceleration at the top and bottom of the defect caused by the viscous diffusion is 
much weaker than that produced at the sides by the lateral divergence of the —u'v' Reynolds 
stress, induced by the pancake-eddy turbulence (cf. figure 6). (The lateral asymmetry of the 
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Figure 12. Terms in mean streamwise momentum budget at Nt = 18 (NEQ regime): (a) rate 
of change du/dt ; ( b ) — duw/dz ; (c) vertical advection — wdu/dz ; (d) lateral advection 
— vdu/dy ; (e) vertical Reynolds-stress divergence — du'w'/dz ; (/) lateral Reynolds-stress 
divergence — du'v'/dy ; (g) vertical viscous-stress divergence v d 2 uj dz 1 . Lateral viscous- 
stress divergence v d 2 u/dy 2 is negligible compared to terms shown. Terms given in units of 
U 2 / L r , where U r = 2.09 Udo and L r — h Q . 
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Figure 13. Terms in mean streamwise momentum budget at Nt = 82 (Q2D regime): (a) rate 
of change du/dt ; (6) — duw/dz ; (c) vertical advection — wdu/dz ; (d) lateral advection 
— vdu/dy ; (e) lateral Reynolds-stress divergence —du'v'/dy; (/) vertical viscous-stress di- 
vergence vd 2 u/dz 2 . Vertical Reynolds-stress divergence —du'w'/dz and lateral viscous- 
stress divergence vd 2 ujdy 2 are negligible compared to terms shown. Terms given in units 
of U 2 /L r , where U r = 2.09 Udo and L r — h Q . 
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Figure 14. Terms in mean streamwise momentum budget at Nt = 586 (late Q2D regime): 
(a) rate of change du/dt ; (6) lateral Reynolds-stress divergence —du'v'/dy; (c) vertical 
viscous-stress divergence isd 2 u/dz 2 . Vertical mean- flow momentum divergence — duw/dz , 
advection —vdu/dy — wdu/dz : vertical Reynolds-stress divergence —du'w'/dz and lateral 
viscous-stress divergence vd 2 ujdy 2 are negligible compared to terms shown. Terms given 
in units of U 2 / L r , where U r = 2.09 Udo and L r — h Q . 


—du'v'/dy contours in figure 145 is a statistical anomaly, due to the relatively small number 
of eddy structures included in the streamwise average at this time.) The decoupling of the 
viscous/vertical and turbulent/lateral mixing mechanisms during the late-Q2D regime is 
consistent with the disparate rates of vertical and lateral spreading observed in figures 7 (c, d ) , 
with /3 Z larger than /3 y for Nt >100, and the earlier observation that the growth of h z is a 
purely viscous phenomenon during this time. 

3.5 Turbulence kinetic energy history and budget 

Figure 15a summarizes the manner in which stratification alters the cross- wake- averaged 
potential and kinetic energy of the turbulence, ^g(p' p') / poo\T\ (symbols) and (k) = 

(lines), respectively. (Recall the cross- wake average is over regions where u > 0.05 Ud •) 

The density fluctuations are largest early in the NEQ regime, when buoyancy most 
strongly affects the flow, and the internal waves are the most active. The present case, 
with its much larger initial Froude number (F^ ~ 130), does not contain buoyancy-induced 
temporal oscillations of TKE at early times. This is in contrast to Dommermuth et aid s 
[8] Foq = 4 LES and Brucker & Sarkar’s [9] = 8 DNS, for which dk/dt first becomes 

positive near Nt = 5.5 for the former and Nt = 2.3 for the latter. (See figure 21 regarding 
values of the virtual origin Nt vo used for those simulations.) Those oscillations, which are 
damped more slowly in the lower F 0 0 flow, are indicative of transfer back and forth between 
the potential and kinetic energy of the turbulence, associated with changes in sign of the 
TKE-to-turbulence potential energy (TPE) exchange term, B t k e = —gp'w'/pQQ (Brucker & 
Sarkar mention, but do not quantify, ‘large values of the buoyancy flux at early time’). For 
the present flow, B t k e always acts as a net sink of TKE (figure 35) despite the monotonic 
decrease with time of the local mean Froude number Fa (figure la, 5). 
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Figure 15. Histories of (a) cross- wake-averaged turbulence potential and kinetic energy, and 
(5) Reynolds-stress anistropy bij = (lu^u'jl) /2(k) — Sij/ 3. Inset plot in a shows components 
scaled by local U a • 


Instead of increasing towards the beginning of the NEQ regime, here the TKE increases 
for 30 < Nt < 70, i.e. during the NEQ-to-Q2D transition defined by the Ua history (fig- 
ure lb). The reasons for this increase, which are unrelated to the TPE-to-TKE transfer via 
H t ke, will soon be clear. 

The difference between the rates of decay of the three velocity fluctuations, as well as the 
non-constant variation of (k)/U% for Nt > 4 (solid curve in figure 15a inset), demonstrate 
the anisotropy introduced by the stratification, and how it breaks the early-time axisym- 
metry, leading to a flow for which the relationship between the magnitude of the velocity 
fluctuations and the maximum mean defect continues to evolve throughout the NEQ and 
Q2D regimes. Of the three components, the streamwise TKE ^(u'u') (dashed curve) ex- 
hibits the least variation with respect to the local Pj, but even so varies significantly (by 
over a factor of five) as the wake regimes come and go. For Nt > 150, when the TKE is 
reasonably approximated by a £ -4 / 3 power law, k decays more slowly than Pj does, such 
that (assuming Ud ~ £ -0 ' 85 i n the Q2D range) k/U\ ~ £ +0 ' 37 (figure 15a inset plot). 

The disparate effect of the stratification upon the three velocity- fluctuation components, 
and thus on the structure of the Reynolds-stress tensor, can also be seen in figure 155. 
These Reynolds-stress-anisotropy histories show that the velocity fluctuations first deviate 
from axisymmetric- decay behavior (i.e. v'v' ~ w'w' and u'v' ~ u'w') near Nt = 2, the time 
at which the mean-defect height and width begin to differ from each other (figure 7c, d). 
Thereafter, the vertical component falls monotonically, first slowly then precipitously, such 
that w'w' and both u'w' and v'w' are negligible for Nt > 40 (cf. Spedding [12] and Diamessis 
et al. [2]). The absence of vertical velocity means that buoyancy can no longer directly cause 
a transfer between the turbulence potential and kinetic energies. 

The streamwise normal stress u'u' becomes a larger and larger fraction of the 2k — 
u'u' + v'v' + w'w' total until the flow is well within the NEQ regime, when dk/dt becomes 
positive, near Nt ~ 30. Afterwards, the normal component v'v' increases faster than u'u ' , 
and u'v' is the only significant shear stress. The dominance of u'v' corresponds to that of 
the Pi 122 = —u'v'du/dy production term; recall the importance of —{Pi 122 } in the {K 11 } 
budget and see below. During the Nt > 70 phase of the Q2D regime, when dk/dt is again 
negative (and the flow is dominated by pancake-eddy structures), the lateral fluctuations 
have more energy than the streamwise fluctuations do, with (v'v') / (u'u') 2 for Nt > 400. 
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Figure 16. History of terms in the budget of integrated TKE {k}. Terms in (a) given in units 
of L r , in ( b ) as fraction of local {e t ke} • Interpolant in a (chain-double-dot curve) given by 
A e (Nt )~ 5 / 3 = A e (Nt— Nt vo )~ 5 / 3 , with ( A e ,Nt vo ) = (0.00138,0.645). Production terms not 
shown are either zero because of homogeneity ({Pun} = {P 2211 } = {P 3311 } = 0) or negligi- 
ble ({P 2233 } ~ {P 3322 } ~ 0). Components of dissipation rate e aa = v(du' a /dxj)(du' a /dxj) 
(no sum on a), where ea = e t k e — i'd 2 u , i u , -/dxidxj , and e t k e = 2 z/sCsC. 


We now consider the transport equation for the TKE k, 
dk 

— A:ke T P tke T F tke T Ptke ^tke T PVtke, (7) 

where the rates of change due to advection Ake, production P t k e , flux divergence Ake, 
buoyancy transport B t ke and dissipation e t k e are, respectively, Ake = — Ujdk/dxj , P t ke = 
Piijj = ~ u'yjdui/dxj , Ake = dfj/dxj , with fj = -\u , j u , i u\ - ^'p'/Poo + vdu'^.jdxi + 
vdk/dxj , and P tke = and e tke = 2 vs' ij s' ij , with sC = \(du'Jdxj Adu^/dxi). 

For the flow at hand, these reduce to Ake = —vdk/dy — wdk/dz , P t ke = Pa 22 + Pzi33, 
Ake = d[— -v'u'iu'i— v'p'/ poQ+vdv'v 1 / dy+i'dv'w , /dz+i'dk/dy\/dy+d[—^w'u' i u' i —w'p'/ P00+ 
vdv'w' / dy -\- isdw'w' /dz + vdk/dz\/dz and P t k e = —gp'w'/p 0 c . The work done by the fringe 
forcing Fi in ( 1 ) is W t k e = u[Fi (i = 1 , 2 , 3 ). 

We begin by examining the cross-plane- area-integrated histories of the terms in (7), 
shown in figure 16. (This information was presented in a different format in figure 35; 
it is repeated here to highlight the relative importance of individual budget terms.) The 
cross-plane averages used in figure 15 can be approximated from the cross-plane integrations 
used in figure 16 by dividing the latter by the cross-sectional area of the mean defect, such 
that e.g. (k) ~ {k}/A7rh y h z . Since the advection and flux- divergence terms integrate to 
zero, figure 16 reveals the contributions of the pure source Ptke, sink e t k e and TKE-TPE 
transfer P t ke- The cross-plane variation of these terms, and Ake(p, z) and Ake(p, z), will be 
investigated below. 

For the present high-Poo flow, the fringe- work contribution {W t k e } to the TKE balance 
is negligible (recall that the production, dissipation and buoyancy sum shown in figure 3b is 
sufficient to provide an excellent representation of the d{k}/dt history, for 0 < Nt < 840.) 
The net buoyancy transfer is also much smaller than the production and dissipation, but is 
not negligible, with —{B t k e } reaching a maximum of about 30% of the integrated dissipation 
{e tke } from early in the NEQ regime (Nt ~ 5) and remaining fixed near this ratio until 
Nt « 70 (figure 165), after which —{B t k e } is a monotonically decreasing fraction of {e t k e } 
until it becomes negligible after Nt « 150. Between Nt ~ 25 and 70, all three non-negligible 
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(c) Nt = 18 (d) Nt = 40 



( y-Vc)/h y ( y-Vc)/h y 


Figure 17. Contours of buoyancy transfer B t k e as a fraction of (e t ke)max, the local maximum 
e t ke from figure 20, at: (a) Nt = 2.5, ( b ) Nt = 5.5, (c) Nt = 18, (d) Nt = 40, (e) 
Nt = 82 and (/) Nt = 586. Minimum (darkest) contour level -Cy im corresponds to the 
most negative B t k e /(ctke)max at each time. White solid-line contours indicate five levels of 
positive #tke/(etke)max, over range +0.01|Ci im | < i3 t k e /(etke)max < [^tke/(etke)max]max, where 
[^tke/(^tke)max]max is maximum ^tke/(^tke)max 5 which at N t — 2.5, 5.5, 18, 40, 82 and 586 
is respectively +0.0007, +0.085, +0.37, +0.14, +0.06 and +0.005. Black solid-line contours 
correspond to u — 0.05 Ud and 0.5 Ud- 
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(a) Nt = 2.5 


(6) Nt = 5.5 




Figure 18. Contours of density- fluctuation variance p' p' at: (a) Nt = 2.5, ( b ) Nt = 5.5, 
(c) Nt = 18, (d) Nt = 40, (e) Nt = 82 and (/) Nt = 586. Values given in units of 
(pooUr / gL r ) 2 . Solid-line contours correspond to u — 0.05 Ud and 0.5 Ud- 
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terms are roughly constant in absolute terms, with d{e t ke}/d£ and d{B t ke}/dt approximately 
zero and d{P t k e }/d£ increasing slightly. Spedding [12] found a similar NEQ plauteau in the 
total dissipation e t k e + Dujj in his high- Poo experiments, which he attributed to ‘vertical 
shear taking over’. We have already observed that the global/integrated buoyancy transfer 
{Btke} nets as a pure sink throughout the lifetime of the present flow (figure 36). At times, 
however, there are local regions in which B t ke is positive, where the potential energy of 
the turbulence is converted to kinetic energy. During the early NEQ phase, these tend 
to be concentrated just above and below the wake (see white-line contours in figure 176), 
while later they are more diffuse and disorganized, although still occurring mostly outside 
the wake core (figure 17c-e). Much later, in the viscous-dominated Q2D regime, when 
{Btke} / {^tke} — > 0, positive Btke tends to lie in the non- vortical regions to the left and right 
sides of the wake (figure 17/). The relationship between B t k e and the density fluctuations 
p' p' can be inferred by comparing figures 17 and 18. The regions of maximum negative B t ke 
are in general highly correlated with those of maximum \p'\, although the departure from 
axisymmetry associated with the beginning of the NEQ regime is more apparent in the p’ p' 
field at Nt = 2.5 (figure 18a) than in the corresponding B t ke distribution (figure 17a). The 
local p’ p' maxima at Nt = 2.5 at the top and bottom of the wake, near z = ±1.5 h z , are 
consistent with the ‘excess’ TPE that later ‘springs back’ into TKE at these locations, by 
way of the positive B t ke observed in figure 17 b at Nt = 5.5. 

Figure 16a shows that the Reynolds number is large enough for the early/3D-regime 
evolution to satisfy the inertial scaling e t k e ~ /c 3 / 2 /h, such that e t k e ~ £ -7 / 3 and {e t ke} ~ 
£ -5 / 3 (chain-double-dot curve; cf. Brucker & Sarkar [9]). For late times, approximately 
Nt > 200, the decay of {Vtke} and especially {etke} is reasonably well approximated by t~ n 
with n = 1.7 ~ 5/3. Thus, paradoxically, the TKE dissipation evolves according to the 
inertial scaling during both the 3D and Q2D regimes, despite their qualitatively different 
eddy structures (see figure 20 d). 

Once the wake turbulence is fully established, the net rate of TKE dissipation is larger 
than that of the production V t ke , except for a substantial period during the NEQ regime, 
from Nt ~ 20 to 90, when {Vtke} is significantly larger than {e t ke} (and gradually increases, 
in absolute terms), after which it returns to its pre-NEQ behavior, with {Vtke} / {^tke} < 1 for 
the Q2D phase. This explains the growth of TKE during this period, mentioned above. From 
figure 166 we find that the enhanced NEQ-regime production is due solely to amplification 
of the Pi 122 = —u'v'du/dy component (dotted curve), since the —u'v' stress grows more 
rapidly than du/dy decreases (and the fall of u'w' is more rapid than the attenuation 
of the decay of du/dz\ see figure 166 and histories of the MKE dissipation components 
P1122 = v(du/dy ) 2 and P1133 = v{du/dz ) 2 in figure 10). This ‘structural’ alteration, 
represented by the decrease of \u'w'\ and increase of \u'v'\, was also found by Brucker & 
Sarkar [9] in their DNS for both towed and self-propelled cases. They concluded that 
stratification causes the TKE production to mainly consist of Pi 122 during this time. The 
present results underline the extent to which this is true, and that in fact P t k e can effectively 
be replaced by Pi 122 for ah times after Nt « 20. 

While both the present and Brucker & Sarkar simulations contain enhanced production, 
only here does it lead to growth of TKE (and not just reduction in the net TKE decay) 
during the NEQ regime. This is because of the large amount of energy radiated by internal 
waves between Nt = 25 and 75 in Brucker & Sarkar’s Poo = 8 flow, which counteracts the 
impact of the Pi 122 growth. The strength of the internal waves is an indication of when 
the TKE growth occurs. When they are strong, at low Poo, buoyancy leads to an early 
TKE increase driven by the TPE-to-TKE transfer, before the NEQ regime begins - during 
which the energy the waves radiate mitigates the effect of the enhanced TKE production 
associated with the flattened, two-layer NEQ structure, and d{k}/dt is small but negative. 
When the waves are weak, at high Poo, they do not attenuate the enhanced P t k e ~ P 1122 , 
leading to net TKE growth during the NEQ regime. 


34 


Also included in figure 165 are histories of the three components of the integrated rate of 
homogeneous dissipation, {e^} = {en} + {e 22 } + {£33}, as a fraction of the full (actual, ther- 
modynamic) rate of TKE dissipation {e t ke} = 2i /{s'-s'-}, where eu = u(du f i /dxj)(du f i /dxj). 
(There is very little difference between the integrated homogeneous and full/inhomogeneous 
dissipation for this flow, with 0.9999965 < {e^}/{e t ke} < 1.0012 for all times considered.) 
Examining e aa for a = 1, 2, 3 allows us to assess the interplay between the rates of produc- 
tion and dissipation of the individual normal-stress components u' a u' a . (Other dissipation 
decompositions have been utilized by Spedding [ 12 ] and Brucker & Sarkar [9]; see also 
Itsweire et al. [33].) 

The vertical component {e 33 } ~ {(dw' / dx) 2 } + {(dw' / dy) 2 } + {(dw' / dz) 2 } falls mono- 
tonically as w' does, with d{e 33 }/d t becoming negative from Nt ~ 7 until Nt > 50, when 
{ 633 } is negligible with respect to {e t ke} • 

For the streamwise component, the increase in Pi 122 for Nt > 5 is accompanied by an 
increase in {eu}, such that the net effect (to which the pressure-strain correlation 0 n = 
p'(du' /dx)/ poo may also contribute) is to arrest the decay of u'u' , such that d{u'u'}/dt first 
changes sign and then becomes slightly positive as the flow enters its NEQ phase (figure 15a). 
This is an example of what is often referred to as an ‘extra strain’ effect (Bradshaw [34]), 
which involves a difficult-to-model change that results from the imbalance of much larger 
terms of opposite sign, each of which change much more than the net imbalance. Similar 
behavior is found in the Q 2 D regime, in that as d{u'u'}/dt becomes negative again, both 
{P 1122 } and {eu} decrease (the latter accounting for roughly one-third of the total {e t ke}) - 

The evolution of the lateral component during the 3D-to-NEQ transition is not directly 
affected by production (since P 2211 = 0 and P 2222 ~ P 2233 ~ 0). Consequently, the slight 
increase in {v'v'} between Nt « 25 and 80 shown i n figure 15 a must be due to inter- 
component transfer via the pressure-strain term </> 22 = p'(dv' /dy)/ poo, mitigated by {£ 22 }, 
which first falls slightly below {etke}/3 (i.e. the isotropic value approximated by all three 
components before Nt ~ 5), then grows to account for about 2/3 of the total. 

The dissipation histories are broadly consistent with those observed in Brucker & Sarkar’s 
[9] low- P qo / self-propelled wake simulation, in that (du'/dz) 2 and (dv' /dz) 2 (i.e. the domi- 
nant terms in eu and 622 , respectively) are respectively the second-largest and largest terms 
for Nt > 10 (although the contributions of dv' /dx, dw' /dx and dw' /dz were not quan- 
tified, these are likely to be negligible). The {e aa } histories are thus also consistent with 
the enhanced vertical gradients associated with the flattened eddy structures of the NEQ 
and Q 2 D regimes (figures 5 and 6 ). We note that {^ 22 } / {^n } ~ 2 for Nt > 150 while 
(v'v')/ (u'u') > 1.7 in the late-Q 2 D regime, i.e. for Nt > 500 (figure 155). 

Terms in the TKE budget for the (essentially) unstratified axisymmetric benchmark at 
Nt = 1.2 are illustrated in figure 19. Shown are radial profiles of the production P t k e , 
dissipation e t k e , buoyancy transport P t ke (which at this time is very small) and the flux 
divergence P t ke, which solely redistributes TKE within (and thus integrates to zero over) 
the y-z plane, and thus played no role in the integrated budget histories considered to 
this point. The importance of the flux divergence (open symbols) is indicated by the large 
contribution it makes to the net dk/dt (solid symbols), especially to the TKE growth at 
the edge of the wake. The shapes of the production and dissipation profiles are the same as 
those found in the self-similar axisymmetric experiment of Uberoi & Freymouth [35], but 
here the transport increases the TKE not just at the edge of the wake (as it does in the 
experiment) but also near the centerline. (The difference is perhaps a consequence of all the 
fluxes, including the viscous ones, being included in P t k e in figure 19, rather than only the 
inviscid, boundary-layer terms shown in [35].) But while the profile shapes of the budget 
terms for the Uberio & Freymuth experiment and the present simulation are comparable, 
their relative magnitudes are not - which provides further evidence for the non-universal 
nature of the turbulence for (unstratified) wakes generated by different objects [5,36],RCC. 

The influence of the stratification upon the three dominant terms in the TKE budget, 
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Figure 19. Radial profiles of terms in the TKE budget at Nt = 1.2: , rate of production 

7^tke5 5 (minus) rate of dissipation — e t ke5 , rate of buoyancy transfer B t k e ; °, flux 

divergence Vtke] ♦, 'Ptke — ^tke + Btke + *Ftke- Advection —vdk/dy — wdk/dz is negligible 
compared to terms shown. Budget terms normalized by U^/L r , where reference velocity 
U r = 2.09 Udo and reference length L r — h Q . 


and in particular how they evolve in the y-z plane, can be observed in figure 20. To facilitate 
comparison, ‘quadrant- averaged’ contours of dk/dt , Vtke , — e/c and Vtke at selected times are 
combined in each plot, in the ‘northwest’, ‘northeast’, ‘southwest’ and ‘southeast’ quadrants, 
respectively. The lateral and vertical coordinates are scaled by their respective width or 
height, to clarify the spatial variations, and reveal the early ‘kinematic’ stratification effects, 
which act solely to flatten and widen the axisymmetric distributions. The black lines trace 
the shape and extent of the mean defect, via 0.05 Ud and 0.5 Ud contours, while the white-line 
contours indicated the dominant term(s) at each time. 

At Nt = 2.5 (figure 20a), the only indication of the beginning of the end of the 3D regime 
is the slight lack of axisymmetry in the flux term, which causes the net TKE growth at the 
edge of the wake to be larger at the sides than at the top and bottom. The production 
and dissipation are not yet deeply affected, since using the width and height to normalise 
the cross-plane axes causes them to maintain their (axisymmetric) radial shapes from the 
unstratified state (figure 19). The location of the maximum production is still just inside the 
u = 0.5 Ud (i.e. r = h) contour (close to the location of the maximum TKE; see figure 85), 
and the maximum dissipation falls between r — h and the centerline. As is the case for 
the unstratified/axisymmetric state, the flux term plays a large role in setting the net 
dk/dt variation, and is solely responsible for the TKE growth at the wake edge (since the 
dissipation tends to cancel the production here). 

When the flow is well into the NEQ regime, tF t ke is much less significant, leaving the 
production-dissipation imbalance to control dk/dt (figure 205, for Nt = 18). The dissipation 
has by this time become quite anisotropic, with e t k e ~ 2en ~ 4e22 ~ 5e33 (figure 165). 
Nevertheless, its cross-plane variation - while tending to align into two flat regions above 
and below the y axis (consistent with the double-layer structure observed in figure 5, for 
Nt ~ 20) - remains more axisymmetric than that of the production. In contrast, Vtke 
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(a) Nt = 2.5 



(c) Nt = 82 



(d) Nt = 586 
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Figure 20. Contours of terms in TKE budget at: (a) Nt = 2.5, ( b ) Nt = 18, (c) Nt = 82 and 
(d) Nt = 586. Streamwise- averaged quantities have been further averaged (‘double folded’) 
across the lateral y and vertical z axes through the wake centerline (y c ,z c ), and plotted in 
separate quadrants, such that the upper left (‘northwest’), upper right (‘northeast’), lower 
left (‘southwest’) and lower right (‘southeast’) region of each plot respectively corresponds 
to the net rate of change dk/dt = A t ke ± Ttke + ^tke + £>tke _ e t ke 5 production T^ke, (minus) 
dissipation — e t ke and flux-gradient transport Jdke- Buoyancy flux term B t k e is shown in 
figure 17. Advection Atke = —vdk/dy — wdk/dz is negligible compared to the other terms. 
Budget terms normalized by U^/L r , where reference velocity U r = 2.09 Udo and length L r = 
h Q . White solid- line contours indicate ±0.67 of the larger of the minimum and maximum 
magnitude at each time. Black solid-line contours correspond to u = 0.05 Ud and 0.5 Ud- 
The (maximum, minimum) values at Nt = 2.5, 18, 82 and 586 are respectively (±1.49 x 
10“ 5 , -2.25 x 10“ 5 ), (±1.77 x lO” 7 , -1.1 x 10” 7 ), (±5.3 x 10” 8 , -3.2 x 10” 8 ) and (±8.15 x 
10 — 11 , —2.6 x IQ -10 ). 
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tends to be concentrated into a ‘wing-like’ region near y = ±h y , because it is now nearly 
all associated with the —u'v' stress (Meunier et al. [3]), with Vtke ~ P 1122 — —u'v'du/dy. 
Consequently, the net dk/dt w Vtke — etke tends to be segregated according to sign, into 
either an ‘hourglass ’-shaped zone aligned with the z axis (TKE decay) or double ‘hotspots’ 
of TKE growth focused near y = ±h y . The latter is compatible with the ‘double-pronged’ 
structure in the TKE contours at Nt « 20, seen in figure 8 c. 

The flow’s entrance into its Q2D phase is characterized by the continuing importance of 
Vtke ~ P\\ 22 i which is again concentrated along the y axis, near y =» ±h y (figure 20 c, for 
Nt = 82). However, the production now contains a small region of negative Vtke centered 
at the origin, due to negative P 2222 — —v'v'dv/dy. This is due to the combination of 
finite v'v' (whose maximum occurs at the centerline), and the positive dv/dy associated 
with the mean stagnation flow (i.e. vertical collapse and lateral expansion) at the origin. 
The dissipation is now even more concentrated in two flat regions above and below the 
y axis, within the u = 0.5 Ud contour, than it was at Nt = 18 - corresponding to the 
enhanced vertical shear mentioned by e.g. Spedding [ 12 ], and the primary roles played at 
this time by qi ~ ( du'/dz ) 2 and 622 ~ (dv' /dz) 2 (figure 166). The flux- divergence term 
has grown in importance, compared to its Nt =18 levels, such that it significantly reduces 
Vtke (near y = Fh y ) and decreases c t ke (near z = ±h z ). The net result leads to the TKE 
distribution shown in figure 8 d, which demonstrates the persistence of the flattened double- 
peaked structure as the flow passes from the NEQ to the Q2D regime (see also Brucker & 
Sarkar [9]). 

Much later, at times for which the turbulence is dominated by the quasi-2D pancake 
eddies displayed in figure 6, the Vtke ~ Pi 122 production is still active near y — ±h y , and e t ke 
(~ 1.5e22) is still focused in horizontal layers near 2 = ±0.5 h z , but both are weaker than the 
flux term Vtke- As is the case for Vtke and e t ke 5 the distribution of Vtke is closely correlated 
with the two-layer pancake-eddy structure at this time. The largest (positive) vertical- 
gradient contribution dfe/dz is aligned with the horizontal layers centered near z = ±h z /2 
(and thus reduces the effect of the dissipation here), while the smallest (negative) values lie 
between the layers, along the y axis; the maximum and minimum lateral divergence df^/dy 
are clustered along the y axis, at the lateral edges of the dual Vtke (and TKE) peaks, with 
the negative values closer to the wake centerline. The flux term combines with Vtke and e t k e 
to yield large negative dk/dt in a shallow horizontal layer near the centerline. While there 
is still some TKE growth due to the combined effect of V t ke and Vtke, between the y and z 
axes along the u = Ud contour, it is very weak. 

3.6 Comparison with previous simulations 

Figure 21 shows the mean-defect histories from earlier simulations at various Froude and 
Reynolds numbers, along with the present results (Fdo = 130, Rdo — 2625). Available 
DNS and a recent LES are included: The DNS of Gourlay et al. [13] ( Fdo = 34, Rdo — 
1472) and Brucker & Sarkar’s [9] Case TR50F04 ( F do « 0.75, R do « 3250) and TR10F20 
(Fdo ~ 3.75, Rdo ~ 650), and Diamessis et aV s [2] Case R100F64 LES ( Fdo ~ 34, nominal 
Rdo ~ 11,000). Also shown are the unstratified/3D-regime power-law diagnostics for each 
case (inset plots in figure 21a, b). These are used to estimate virtual-time origins for each 
simulation, to facilitate comparison by plotting histories in terms of Nt = Nt — Nt vo . 

The various histories cannot be expected to agree in the scaling used for figure 21, which 
takes no account of Froude and/or Reynolds numbers, since the maximum velocities Ud, 
heights t z , and widths t y are presented solely in terms of their initial values. (Also, t z is 
for some cases given in terms of the 0. 5 Ud half-height h z and for others in terms of integral 
measures; see figure 21c caption.) Nevertheless, several common features can be observed, 
including similar power-law behavior over the same Nt range. The height decreases with 
time for Nt > 2 in every case, although the reduction of h z is larger than that experienced 
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Nt Nt 


Figure 21. Histories of (a) maximum velocity ( b ) height-to- width ratio £ z /£ y , (c) height 

£ z and ( d ) width ^ of mean defect: , present; o, [13]; ■, Case TR50F04 of [9]; •, 

Case TR10F20 of [9]; o, Case R100F64 LES of [2] . Inset plots in (a) and (6) show unstratified 
similarity diagnostics for (a) velocity x — ( Ud/Udo)~ 3 ^ 2 and (6) height £ z = (£ z /£ 0 ) 3 (lower 
curve/symbols) and width £ y = (£ y /£ 0 ) 3 (upper curve/symbols), which are used to define the 
virtual-time origin t vo for each case, with Nt vo = 0.645 for the present DNS, Nt vo = 0.2 for 
Brucker &; Sarkar’s TR50F04, Nt vo = —0.5 for Diamessis et al. and Nt vo = 0 for Gourlay et 
al. [13] and Brucker & Sarkar’s [9] TR10F20. For the present DNS, the linear curve-fit for x 

and £ ( in upper right-hand-side inset plots in a and b) are x — ( Ud/Udo )~ 3 ^ 2 = B • Nt 

and £ = (r*/rj) 3 s= C • Nt, where Nt = Nt — Nt vo and {B,C, Nt vo ) = (3.1,2.9,0.645). 
The maximum defect velocity Ud in (a) for Brucker & Sarkar’s [9] Case TR50F04 (but not 
TR10F20) is replaced by their Integral measures are used for the height and width 

for the present DNS, with t z — r*, £ y = r*, and t Q the value of r* z — r * at Nt = 0. For 
Brucker & Sarkar’s [9] Case TR50F04, £ z and £ y are respectively represented by integral 
scales (in their notation) Rs and R 2 , with £ 0 the half-radius of the initial Gaussain mean- 
velocity profile, while for Diamessis et al. [2], and effectively for Gourlay et al. [13], £ z = h z , 
£ y = h y and £ 0 = h Q , where h z and h y are determined from fitting an assumed Gaussian 
profile. 
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by the integral-measure heights (which is consistent with the behavior found for the present 
DNS, in that both dh z /dt and drj/dt are negative for 8 < Nt < 18, but h z falls faster than 
r*; figure 7b). 

The 3D-to-NEQ transition also tends to occur at about the same time for each of the 
simulations, with Nt between 2 to 2.5 the point at which the heights and widths first begin 
to differ from each other; this can be seen in figure 215, in both the height-to- width ratio 
£ z /£ y and is even more obvious in the unstratified/axisymmetric height /width diagnostic £ = 
(£/I 0 f shown in the inset plots (recall that linear variation with time of £ is characteristic 
of the pure-axisymmetric 3D regime). 

The reasonably good collapse of the mean aspect ratio t z j i y over the early NEQ periods 
of the five simulations was not anticipated, both because of the wide Froude and Reynolds 
number ranges, as well as the differences in the initialization strategies used for the various 
simulations (e.g. Brucker & Sarkar’s [9] DNS held the mean velocity fixed until Nt = 
0.2 (TR10F20) and 1 (TR50F04), while Diamessis et aids [2] LES gradually introduced 
the background stratification during a precursor simulation). The £ z /£ y histories again 
demonstrate the increase of the Nt range of the NEQ regime with Reynolds number - 
compare the times at which the £ z /£ y histories begin to deviate from the negative power- 
law behavior for NEQ, for Brucker & Sarkar’s [9] TR10F20 (40 < Nt < 100), Gourlay et 
al. [13] (Nt ~ 70), present (Nt ~ 105), and Brucker & Sarkar’s [9] TR50F04 (Nt > 400, at 
least). (The late-time I z /I y history for Diamessis et al. [2] (open diamonds) is thought to be 
statistically unreliable, in light of the relatively small streamwise domain used.) The same 
message regarding the Reynolds-number dependence of the NEQ period can be inferred 
from figure 21c, by examining the durations of the dl z /dt « 0 regions for each case. 

The tendency for i z and £ y to deviate from the 3D/axisymmetric ideal sooner than 
Ud does, which was seen above (figure 7), also holds for the other simulations (although 
not as markedly), since the unstratified/axisymmetric similarity diagnostic for Ud, X — 
(Ud/Udo)~ 3 / 2 , remains linear longer (until Nt ~ 2 to 2.5) than the corresponding height/width 
diagnostic, £ (Nt w 2 to 5); cf. figures 21a and b inset plots. 

A scaling to account for the significant Froude-number sensitivity of the Ud histories 
seen in figure 21a has been developed by Meunier & Spedding [17]. It amounts to assuming 
Ud will eventually depend solely on the volume- flux deficit Id, the Brunt- Vaisala frequency 
N and time t, which leads via dimensional analysis to 

- Nto), ( 8 ) 

where 4> is a universal non-dimensional function of Nt — Nto and Nto is a non-universal 
offset time (unrelated to the 3D virtual-origin time Nt vo ), which will depend for example 
on Reynolds number and the near- field/early-time wake evolution. In terms of a virtual 
wake-creating sphere of diameter D emersed in a steady uniform stream of density p ^ 
and speed Vqo, experiencing a drag force Fd, we can write Id = \V oq D 2 Cd, where Cd = 
Fd / (K pooV^oD 2 /&). Equation (8) then implies 

= f(C D )*(Nt - Nto), (9) 

^OO 

where /(Cd) — (kCd /2) 1 / 3 . Equation (9) shows the limitation of attempting to collapse 
Ud/Voo data from wakes generated by various bluff bodies (with variable Cd) by using only 
F 0 o (e.g. by plotting (Ud/V^F^ 3 versus Nt', cf. [1,9]). Following Meunier & Spedding [17], 
(8) can be recast in terms of the virtual effective diameter D e q = D(Cd/ 2) 1 / 2 for bluff bodies 
of arbitrary shape, 

TT ¥ l« = ” 1/3 HNt - Nto), (10) 

^OO 
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(a) 





Figure 22. Histories of rescaled mean defect (a,c) as Q2D power- law/ similarity diagnostic 
and (6, d) with respect to Q2D time shift Nto, assuming (a, b) n = 3/4 and (c, d) n = 0.85: 

, present DNS; o, Gourlay et al. [13]; ■, Case TR50F04 of Brucker & Sarkar [9]; 

• , Case TR10F20 of [9]; o, Case R100F64 LES of Diamessis et al. [2]. Thin-solid curves 

( ) in (a, c) are linear curve fits of Q2D diagnostic (which for [13] and Case TR10F20 

of [9] are nearly coincident), used to define time shift Nto in (5, d), such that for ( b ) (n = 3/4) 
Nt 0 =: 19 (present DNS), -55 ( [13]), -63.5 (Case TR50F04 of [9]), -57 (Case TR10F20 
of [9]) & 11 (Case R100F64 of [2]), and for ( b ) (n = 3/4) Nt 0 =: 7 (present DNS), -70 
( [13]), -66.5 (Case TR50F04 of [9]), -71 (Case TR10F20 of [9]) & -3 (Case R100F64 
of [21). 
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where F e ff = ZVoq/NDqr is the effective Froude number. 

The validity of (8) is tested in figure 22. Meunier & Spedding [17] found the universality 
under consideration cannot be expected until the flow reaches the Q2D regime (see their 
figure 5 b). Consequently, we assume <h(AT — Nto) = Cq2d(AT — Nto)~ n , where Cq2d is 
a universal constant and n ~ 3/4. Similar slopes of the (L^/T^A/' 2 / 3 ) -1 / 7 ' 1 diagnostic are 
found for all five simulations, for both n = 3/4 (figure 22a, with Cq 2 d = 0.08 -3 / 4 « 6.65) 
and the n = 0.85 suggested by the present DNS (figure 22c; Cq 2 d = 0.055 -0 85 ^ 11.8). 
Note the dependence of the beginning of the region of common slope (straight thin-solid 
lines in figure 22a, c), and the offset times Nto (AT-axis intercept), on the Reynolds number 
and initialization scheme used by each simulation. The reasonable collapse of the results in 
figures 22 b and d, produced by using the Nto defined for each case by the corresponding 
curve fits in figures 22 a and c, points to the validity of (8), and thus to that of the founding 
assumptions used in the Meunier & Spedding [17] analysis. 


4 Summary and closing comments 

The DNS results presented above reveal the effect of a weakly stratified, buoyantly stable 
background upon the structure and statistics of a turbulent wake, during its various phases 
of development, including its 3D, NEQ and Q2D regimes. This has been documented by 
the evolution of the mean momentum and energy, turbulence kinetic energy, and the terms 
in their budgets. The dependence of the quality of the statistics on adequate ‘eddy sample 
size’ has been observed, as has the Reynolds-number dependence of the duration, in terms 
of AT, of the NEQ regime. 

The mean momentum budget clarifies the relationship between the turbulence suppres- 
sion due to buoyancy and the turbulent and viscous mechanisms responsible for the vertical 
and lateral growth of the wake. The TKE budget, on the other hand, quantifies the man- 
ner in which the changes in turbulence structure are reflected in changes to the relative 
importance and composition of the terms whose imbalance defines the net rate of change 
of TKE. Surprisingly, despite these changes in structure, the TKE dissipation decays at 
approximately the same rate during the early/3D and late/Q2D regimes. This affirms yet 
again the dissimilarity between the Q2D state and pure 2D turbulence [11]. 

The present results illustrate that the wake collapse, signaling the beginning of the 
NEQ regime, need not be accompanied by the transfer from the potential energy of the 
turbulence to the kinetic energy of the turbulence, as has been proposed elsewhere. Instead, 
the enhanced TKE exhibited during the NEQ regime is solely the result of buoyancy-induced 
structural changes, which lead to increased TKE production. Thus, while the potential-to- 
kinetic energy transfer is observed at early times for lower cases (Dommermuth et al. [8] ; 
Brucker &; Sarkar [9]), it is not an essential feature of the 3D to NEQ transition. 

For the present high-Ebo case, internal wave activity is most pronounced at the start 
of the NEQ regime. However, even during this time, internal waves make a negligible 
contribution to the turbulence energy budget. This is in contrast to the strongly stratified 
cases studied elsewhere, for which the amount of turbulence energy radiated by internal 
waves during the NEQ regime can be comparable to that produced by the mean shear or 
lost via viscous dissipation (cf. Brucker & Sarkar [9]). 

The role of stable stratification in driving the wake towards a universal state has also 
been displayed, in that the Q2D structure is independent of the details of the wake initializa- 
tion (a fact frequently observed in other experiments and simulations). The present study 
also indicates that during the Q2D period, the history of the maximum mean defect Ud (ap- 
propriately normalized, in terms of the initial volume- flux deficit Id and the Brunt- Vaisala 
frequency AT) is a unique function of nondimensional time Nt (relative to a flow- dependent 
virtual offset/origin) . This supports the validity of the framework proposed by Meunier & 
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Spedding [17] in their universal model of the mean-defect evolution for a stably stratified 
towed wake. 
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